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SECTION  I 
INTRODUCTION 


1.  BACKGROUND 

The  Air  Force  Armament  Laboratory  (AFATL)  is  developing  autonomous  seekers  to 
attack  both  mobile  land  targets  and  fixed,  high-value  land  targets.  The  purpose  of  the 
mathematical  structure  described  in  this  report  is  to  aid  the  efficient  development  of  such 
seekers.  This  mathematical  structure — known  as  the  AF ATL  Image  Algebra — has  been 
specifically  designed  for  the  concise  expression  and  clear  representation  of  image  processing 
and  pattern  recognition  techniques.  This  structure  provides  a  common  mathematical 
environment  for  target  detection,  algorithm  development,  optimization,  comparison,  coding, 
and  performance  evaluation.  In  addition,  the  Image  Algebra  provides  a  mathematical  basis 
for  a  universal  image  processing  language  which,  when  properly  implemented,  will  greatly 
increase  a  researcher’s  productivity  as  programming  tasks  required  to  compute  image 
transformations  are  greatly  simplified  due  to  the  replacement  of  large  blocks  of  code  by  con¬ 
cise  algebraic  expressions. 

Several  previous  attempts  to  develop  a  unified  algebraic  approach  to  image  processing 
have  met  only  partial  success  in  expressing  all  transformations  of  gray  value  images,  Refer¬ 
ence  1.  In  contrast,  in  addition  to  meeting  the  design  specifications  mentioned  in  the  previ¬ 
ous  paragraph,  the  AFATL  Image  Algebra  provides  a  complete  unified  algebraic  structure 
capable  of  expressing  all  image-to-image  transformations.  The  Image  .\lgebra’s  foundation 
evolved  from  a  33-month  Air  Force/Defense  Advanced  Research  Project  Agency  (DARPA) 
sponsored  research  effort  known  as  the  Image  Algebra  Project. 

2.  ORGANIZATION  OF  DOCUMENT 

This  report  presents  the  accomplishments  of  the  Image  Algebra  Project  and  provides  a 
tutorial  overview  of  the  AFATL  Image  Algebra.  Section  I  provides  the  background  to  the 
project.  Section  II  provides  a  synopsis  of  the  accomplishments  of  Phase  I  and  Phase  II  by 
tasks.  Sections  III  and  IV  describe  the  mathematical  structure  of  the  Image  Algebra  in  some 
detail.  Section  V,  the  final  section  of  this  report,  discusses  Image  Algebra  software  develop¬ 
ment. 


1 


SECTION  n 

PROJECT  ACHIEVEMENTS 

1.  PHASE  I  ACCOMPLISHMENTS  BY  TASKS 

The  Image  Algebra  Project  was  divided  into  a  two-phase  effort  with  the  first  phase 
spanning  12  months  and  the  second  21  months.  The  main  effort  of  the  first  phase  focused  on 
identifying  the  basic  operands  and  operators  of  the  algebra  and  to  provide  a  first  rough  cut  of 
the  overall  mathematical  structure  defining  the  algebra.  Specifically,  the  tasks  and  achieve¬ 
ments  of  the  first  phase  can  be  summarized  as  follows. 

a.  Consolidation,  Classification  and  Description  of  Techniques 

The  task  of  consolidating,  classifying,  and  describing  existing  image  processing 
transforms,  measurement  techniques,  and  feature  analysis  techniques  was  completed  and  the 
results  were  submitted  in  the  form  of  the  203-page  Image  Algebra  technical  report,  Refer¬ 
ence  2.  This  document  collects  in  one  place  descriptions  of  more  than  100  existing  image 
processing  procedures.  F rom  this  collection,  frequently  occurring  elemental  operations  were 
identified  for  use  in  accomplishing  the  next  task. 

b.  Identification  and  Definition  of  Elemental  Operations 

The  question  as  to  which  operations  and  operands  ought  to  be  regarded  as  elemen¬ 
tal  was  addressed  throughout  the  Phase  I  effort.  Except  for  some  slight  refinements  in 
definition  and  notation,  the  Phase  I  set  of  elemental  operations  and  operands  underwent  no 
major  changes  during  Phase  II  of  the  project.  The  definitions  of  operands  and  operations  can 
be  found  in  Section  III  of  this  report.  Potential  structures  which  include  these  operations 
were  investigated  for  use  in  accomplishing  the  next  task. 

c.  Identification,  Description,  and  Evaluation  of  Mathematical  Structures 

Criteria  for  a  set  of  operations  and  criteria  for  a  mathematical  structure  were  com¬ 
pleted  and  reported  in  Reference  3.  These  criteria  can  also  be  found  in  the  Phase  I  Final 
Report. 

The  identification,  demonstration,  and  evaluation  of  image  .\lgebra  structures  and 
properties  which  are  based  on  the  operations  and  structures  already  identified  were  initially 
reported  upon  at  the  first  program  review  meeting.  The  Phase  I  Final  Report  also  includes 
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this  work.  In  addition,  it  evaluates  a  number  of  mathematical  structures  which  support  the 
Image  Algebra  properties. 

The  demonstration  of  the  capability  and  versatility  of  the  optimal  Image  Algebra 
structure  was  completed  and  reported  in  Reference  4.  Translations  of  numerous  image  pro¬ 
cessing  techniques  into  the  Image  Algebra  and  their  expression  in  terms  of  the  defined  ele¬ 
mental  operations  are  incorporated  in  the  final  report  of  Phase  I. 

2.  PHASE  n  ACCOMPLISHMENTS  BY  TASKS 

The  major  effort  of  Phase  II  focused  on  fine  tuning  the  operands,  operators,  and  the 
oveiall  algebraic  structure  defined  during  Phase  I.  Development  and  implementation  of 
Image  Algebra  software,  and  sponsor  requested  documentation  were  additional  major 
required  tasks.  The  specific  Phase  11  tasks  and  their  accomplishments  are  summarized 
below.  The  first  three  tasks  are  concerned  with  the  mathematical  development  of  the  Image 
Algebra.  They  are  closely  interrelated  and — from  both  a  practical  and  theoretical  point  of 
view — form  a  coherent  unit.  For  this  reason  we  shall  discuss  their  accomplishments  collec¬ 
tively  and  view  them  as  the  accomplishment  of  one  major  task.  Tasks  numbers  coincide 
with  the  sponsor’s  statement  of  work  (SOW)  labeling  scheme. 

a.  Task  4.2.1:  Extend  Properties  and  Relationships 

b.  Task  4.2.2:  Identify  Principal  Properties  and  Relationships 

c.  Task  4.2.3:  Consolidate  Theorems  and  Proofs. 

Our  efforts  in  accomplishing  these  three  tasks  were  extremely  successful.  One  doc¬ 
toral  dissertation,  three  masters  thesis,  and  over  20  conference  and  journal  articles  have 
resulted  from  these  efforts. 

Although  many  basic  properties,  relationships  and  theorems  were  established  during 
Phase  I  of  the  project,  the  main  body  of  the  final  structure,  consisting  of  a  multitude  of  new 
definitions,  relationships,  theorems,  and  pertinent  new  insights  was  completed  during  Phase 
n.  For  example,  the  extension  of  the  template  definition  to  include  parameterized  templates 
and  the  extension  of  the  single  valued  image  algebra  established  during  Phase  1  to  a 
multidata/multiscnsor  image  algebra  were  all  accomplished  during  Phase  II  of  the  project. 
The  final  derived  algebraic  structure  not  only  meets  all  the  sponsor  requested  requirements 
and  specifications,  but  in  many  instances  surpasses  these  requirements  and  specifications. 
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Sections  m  and  IV  of  this  report  describe  the  structure  in  full  detail,  listing  all  principal  pro¬ 
perties,  relationships,  theorems,  various  applications  to  image  processing,  and  the  formula¬ 
tion  of  image  processing  techniques  into  the  language  of  this  structure. 

d.  Task  4.2.4:  Summarize  the  Structure’s  Advantages  and  Disadvantages 

We  first  summarize  the  structure’s  advantages.  The  basic  advantageous  properties 
of  the  algebra  have  been  established  in  the  professional  literature  and  can  also  be  ascertained 
from  the  structure’s  description  in  Sections  III  and  IV.  These  properties  are  best  summarized 
as  follows: 

(1)  its  basic  operands  are  capable  of  modeling  any  image; 

(2)  its  elemental  operations  can  be  combined  to  express  any  gray-level  image 
transformation; 

(3)  its  elemental  operations  are  small  in  number,  translucent  and  simple,  and  easily 
taught  to  potential  users; 

(4)  its  theorems  enable  the  simplification  and  optimization  of  algorithms  by  means 
of  identities  involving  the  operations; 

(5)  its  notation  provides  a  deeper  understanding  of  image  manipulation  operations 
and  is  capable  of  suggesting  new  techniques; 

(6)  its  notational  adaptability  to  programming  languages  allows  the  substitution  of 
extremely  short  and  concise  Image  Algebra  expressions  for  equivalent  blocks  of  code,  and 
therefore  increases  programmer  productivity; 

(7)  the  (standardized)  notation  allows  the  use  of  libraries  of  transformations,  which 
in  turn  lower  algorithm  development  cost; 

(8)  the  algebraic  structure  is  machine  independent  and  language  independent; 

(9)  it  is  applicable  to  image  processing  as  implemented  on  any  machine,  whether  a 
massively  parallel  processor  or  a  sequential  processor;  and 

(10)  the  set  of  algebraic  operations  is  naturally  coherent.  Groups  of  operations  are 
interrelated  by  algebraic  identities  and  inequalities.  There  are  notions  of  associativity,  com¬ 
mutativity,  distributivity,  inverses,  identities,  etc. 
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A  highly  dissonant  (as  opposed  to  coherent)  set  of  operations  would  have  required  more  elar 
borate  notation  for  expressing  algorithms,  resulting  in  a  loss  of  translucency  and  user- 
friendliness.  The  more  naturally  coherent  the  set  of  operations,  the  easier  it  is  to  find  suc¬ 
cinct  notation  for  expressing  these  operations.  In  conjunction  with  Properties  (1)  and  (2),  it 
is  important  to  note  that  we  have  been  able  to  formally  prove  that  this  set  of  operations  is 
sufficient  for  expressing  all  computable  image-to-image  transformations. 

In  regard  to  Property  (6),  a  large  variety  of  well  known  image  processing  tech¬ 
niques,  taken  from  standard  textbooks  and  journals,  has  been  translated  into  the  Image  Alge¬ 
bra.  Some  examples  of  these  techniques  and  their  translations  are  given  in  this  document 
while  a  more  complete  listing  can  be  found  in  the  final  report  of  Phase  I  of  this  project.  The 
usual  formulations  of  these  techniques  are  in  terms  of  word  description,  illustration  and/or 
mathematical  formulae.  In  practice,  however,  a  programmer  would  have  to  reformulate  these 
descriptions  and  equations  into  a  computer  program  which  would  result  in  far  more  lines  of 
computer  code  than  lines  of  description.  Thus,  the  standard  formulations,  when  translated 
into  a  higher  level  computer  language,  are  actually  far  more  involved  than  they  appear.  The 
program  codes  representing  the  simple  mathematical  equation  of  the  discrete  Fourier 
transform  or  Sobel  edge  detector  are  prime  examples  of  this.  In  contrast,  the  translated  ver¬ 
sions  presented  in  this  document  are,  for  all  practical  purposes,  coded  programs. 

Image  Algebra  operations  have  been  implemented  in  the  form  of  a  preprocessor  for 
FORTRAN.  The  input  to  the  preprocessor  is  a  FORTRAN  program  that  includes  the  Image 
Algebra  operands  and  operations.  The  preprocessor  reads  such  a  program  and  converts  each 
Image  Algebra  statement  into  a  sequence  of  FORTRAN  code  that  implements  the  Image 
Algebra  statement.  The  brevity  and  translucency  of  these  programs  imply  the  attainment  of 
at  least  two  of  the  above  mentioned  properties,  namely  the  increase  of  a  programmer’s  pro¬ 
ductivity  due  to  simplification  and  reduction  of  code,  and  the  transparency  of  the 
programmer’s  code  since  each  algebraic  statement  depicts  the  computation  on  the  pixel  level. 
For  example,  a  Government  supplied  Autonomous  Target  Recognition  (ATR)  algorithm  for 
tank  detection  in  Forward  Looking  Infra  Red  (FLIR)  images  translated  into  the  Image  Alge¬ 
bra  and  implemented  via  the  preprocessor,  resulted  in  over  75  percent  reduction  of  code  for 
some  of  its  subroutin  js.  Moreover,  a  significant  number  of  undergraduate  computer  science 
and  electrical  engineering  majors,  assigned  senior  year  research  projects  concerned  with 
developing  various  image  processing  algorithms,  provided  a  convincing  argument  for  Pro¬ 
perty  (3).  Without  any  previous  image  processing  experience,  these  students  were 
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performing  all  the  usual  image  processing  tasks  such  as  image  enhancement,  edge  detection, 
image  coding,  region  labeling,  size  discrimination,  etc.,  within  the  first  2  weeks  into  their  pro¬ 
jects.  These  were  not  canned  programs,  but  written  from  scratch  in  terms  of  Image  Algebra 
notation. 

Because  of  its  simplicity  and  translucency,  one  may  doubt  that  the  Image  Algebra 
can  provide  the  necessary  depth  and  breadth  that  may  be  required  by  future  technology. 
However,  such  doubts  are  easily  dispelled  once  it  is  realized  that  one  subalgebra  of  the  full 
Image  Algebra  is  isomorphic  (mathematically  the  same)  as  linear  algebra,  the  most  used  alge¬ 
bra  in  all  the  sciences.  Another  subalgebra  of  the  Im^e  Algebra  includes  the  algebra  of 
mathematical  morphology,  an  area  of  image  processing  that  has  led  to  many  advances  in 
computer  architectures  and  novel  image  processing  techniques.  Similarly,  some  other 
subalgebras  of  the  Image  Algebra  have  already  yielded  several  novel  image  processing  tech¬ 
niques.  In  particular,  we  have  demonstrated  that  the  Image  Algebra  is  useful  as  a  model  of 
parallel  image  processing  and  as  a  tool  for  the  development  of  parallel  algorithms  for  com¬ 
puting  linear  image  to  image  transforms.  We  have  shown  how  the  algebraic  relationships 
resulting  from  the  structure  of  the  Image  Algebra  can  provide  various  useful  results  and  tech¬ 
niques.  Thus,  for  instance,  we  have  shown  how  the  Image  Algebra  provides  an  alternative 
algebraic  formulation  of  linear  image  processing  that  reflects  both,  contemporary  as  well  as 
future  computing  environments.  References  5  and  6.  Specifically,  the  Image  Algebra  can 
serve  as  an  algebraic  model  for  linear  computations  on  computer  networks  based  on  group 
structures,  Cayley  networks,  and  provide  a  link  between  these  networks  and  existing  alge¬ 
braic  structures.  In  short,  many  deep  theoretical  results  and  their  applications  to  image  pro¬ 
cessing  technology  have  already  been  obtained.  However,  in  view  of  the  completeness  proof 
of  the  Image  Algebra,  these  and  future  results  should  not  be  all  that  surprising. 

As  to  the  structure’s  deficiencies,  it  should  be  apparent  from  the  above  discussion 
and  also  from  Sections  III  and  FV  that  the  structure  is  optimal  for  expressing  image-to-image 
and  image-to-real  or  complex  value  transformations.  However,  image-to-symbolic  domain 
transformations  and  symbolic  domain  representation  and  manipulation  as  used  by  image 
understanding  developers  have  not  been  addressed  in  the  Image  Algebra.  Recently,  notation- 
ally  consistent  operations  and  methods  have  been  developed  within  the  Image  Algebra  that 
can  extract  symbolic  representations  from  images.  The  Image  Algebra  scheme  for  chain  code 
extraction  is  a  good  example  of  such  a  method.  In  view  of  this  capability,  it  would  be  desir¬ 
able  to  develop  notation  that  naturally  unifies  the  image  and  symbolic  domains.  An 
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extension  of  the  Image  Algebra  to  the  symbolic  domain  would  establish  a  comprehensive 
notational  standard  for  all  computer  vision  tasks. 

e.  Task  4.2.5:  Study  the  Feasibility  of  using  Artificial  Intelligence  Techniques 

In  our  feasibility  study  we  found  that  the  use  of  artificial  intelligence  techniques  in 
conjunction  with  the  Image  Algebra  can  be  both  internal  or  external.  Internally,  for  exam¬ 
ple,  special  templates  or  special  template  operations  can  be  defined  through  the  use  of 
artificial  intelligence  techniques.  Externally,  intelligence  techniques  can  be  applied  to  an 
analysis  of  algorithms  which  are  expressed  in  the  algebra,  for  example,  the  optimization  of  a 
specific  algorithm,  the  comparative  evaluation  of  two  algorithms,  or  the  development  of  a 
new  algorithm.  Our  findings  were  reported  to  the  sponsor  in  Reference  7. 

f.  Task  4.2.6:  Demonstration  of  the  Algebra’s  Capabilities 

The  Algebra’s  capabilities  have  already  been  explained  and  listed  in  the  previous 
subsections.  F urther  verification  of  its  remarkable  capabilities  can  be  ascertained  from  Sec¬ 
tion  in.  Actual  demonstrations  of  its  capabilities  have  been  given  at  program  review  meet¬ 
ings  and  several  workshops.  Reference  8.  In  addition,  a  Government  furnished  ATR  algo¬ 
rithm  was  translated  into  Image  Algebra  Fortran  resulting  75  percent  reduction  of  code  and 
vast  improvement  in  the  understanding  of  coded  expressions  due  to  the  translucency  of 
Image  Algebra  expressions. 

g.  Task  4.2.7:  Description/Justification  of  the  Image  Algebra  as  a  DOD  Standard 

The  sponsor  has  been  provided  with  the  documentation  and  justification  as  to  why 
the  Image  Algebra  should  serve  as  a  DOD  image  processing  standard,  Reference  9.  Here  we 
briefly  reiterate  two  main  reasons  for  justifying  the  Algebra  as  a  DOD  standard.  In  simplest 
terms,  the  Image  Algebra  provides  for  a  language  which,  if  properly  implemented  as  a  stan¬ 
dard  image  processing  environment,  will  greatly  reduce  research  and  development  costs. 
Government  savings  due  to  commonality  of  language  and  increased  productivity  will  dwarf 
any  reasonable  initial  investment  for  developing  the  Image  Algebra  into  a  standard  environ¬ 
ment  for  image  processing.  In  addition,  since  the  foundation  of  this  language  is  purely 
mathematical  and  independent  upon  any  future  computer  architecture  or  language,  the 
longevity  of  an  Image  Algebra  Standard  is  assured. 

Although  commonality  of  language  and  cost  savings  are  two  major  reasons  for  con¬ 
sidering  the  Image  Algebra  as  a  standard  language  for  image  processing,  there  exist  a 
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multitude  of  other  reasons  for  desiring  the  broad  acceptance  of  the  Image  Algebra  as  a  com¬ 
ponent  of  all  image  processing  development  systems.  Premier  among  these  is  the  predictable 
influence  of  an  Ima^e  Algebra  Standard  on  future  image  processing  technology.  In  this,  it 
can  be  compared  to  the  influence  on  scientific  reasoning  and  the  advancement  of  science  due 
to  the  replacement  of  the  myriad  of  different  number  systems  (i.e.  Roman,  Syrian,  Hebrew, 
Egyptian,  Chinese,  etc.)  by  the  now  common  Indo-Arabic  notation. 

h.  Program  Management  and  Contractor  Recommended  Tasks 

The  remaining  tasks  (SOW  4.3  -  5.0)  were  concerned  with  program  management 
such  as  review  meeting  responsibilities,  cost/schedule  management,  reports,  data  and  other 
contract  required  deliverables.  All  these  were  carried  out  to  the  sponsors  satisfaction.  In 
addition,  a  contractor  recommended  and  sponsor  approved  task  concerned  with  Image  Alge¬ 
bra  software  development  was  initiated.  Section  V  describes  the  accomplishments  and 
results  of  this  task. 


SECTI'  r,  7" 

IMAGE  ALGEBRA  OPERANDS  AND  OPERATORS 

1.  INTRODUCTION 

Throughout  this  document  we  adopt  commonly  accepted  and,  whenever  possible,  stand¬ 
ardized  mathematical  notation.  In  particular,  Z,  R,  C,  amd  denote  the  set  of  integers, 
real  numbers,  complex  numbers,  and  binary  numbers  of  fixed  length  k,  respectively.  The  set 
theoretic  notions  of  empty  set,  element  of,  not  an  element  of,  union,  intersection,  and  subset 
will  be  denoted  by  0,  €,  ^  U,  D,  and  C,  respectively.  We  assume  that  the  reader  of  this 
document  is  familiar  with  the  basic  notions  of  set  theory.  However,  for  sake  of  uniformity  of 
notation  and  completeness  we  briefly  state  several  basic  mathematical  concepts  and  their 
corresponding  expressions  that  are  used  throughout  the  technical  part  of  this  document. 

Intuitively,  a  function  /from  a  set  X  into  a  set  Y,  written  /:X  — ►  Y,  is  a  rule  which 
assigns  to  each  x  €  X  some  element  y  of  Y,  where  the  assignment  of  x  to  y  by  the  rule  /  is 
denoted  by  y  =/(x).  Throughout  much  of  this  document  we  shall  specify  functions  by 
assignment  rules  and  call  the  set  {  (x,y) :  y  =  /(x)  and  x€X  },  where  /:X  —*■  Y,  the  graph  of 
f.  We  will  also  periodically  refer  to  certain  special  properties  and. types  of  functions.  In  par¬ 
ticular,  it  will  be  important  to  distinguish  between  the  following  types  of  functions; 

A  function  /:X  — *  Y  is  said  to  be  one-to-one  (or  1-1)  if  distinct  elements  in  X  have 
distinct  images;  i.e.  if  /(x)  =  /(z)  ,then  x  =*z. 

A  function  /:X  — ♦  Y  is  said  to  be  onto  (or  /  is  a  function  from  X  onto  Y)  if  every 
y€Y  is  the  image  of  some  x€X;  i.e.  if  y€Y,  then 3  x€X  such  that  y  =/(x). 

If  /:X— *^Y  is  one-to-one  and  onto,  then  there  exists  a  function  Y— >  X,  called  the 
inverse  of  f  and  denoted  by  f~^,  with  the  property  that  /“‘(y(x))=xVx€X. 

A  function  /'.X— Y  assigning  all  x  €  X  to  the  same  single  element  y  €  Y  is  called  a 
constant  function. 


Given  a  function  /;X  — ►  Y  and  a  subset  A  C  X,  then  the  function  /  considered  only 
on  A  (i.e.  the  function  {  (x,/(x)) ;  x  €  A  }  is  called  the  restriction  of  /  to  A  and  is 
denoted  by  Thus,  yly^  =  /n(A  X  Y). 

In  the  reverse  direction,  if  A  C  X  and  /:A  — ►  Y,  then  any  function  /:X  — ►  Y  with  the 
property  /l^  =/  is  called  an  extension  of  /over  X  relative  to  Y,  which  will  also  be 
denoted  by  /l^,  .i.e.,  7  =  /l*. 

Two  important  functions  whose  arguments  are  sets  are  the  choice  function,  denoted 
by  the  word  choice,  and  the  cardinality  function  card.  Given  a  set  X  then  choice(X) 
returns  an  arbitrary  element  of  X,  i.e.  choice(X)=xGX.  The  function  card  on  the 
other  hand  counts  the  number  of  elements  in  X,  i.e.  card(X)=n  where  n  is  either  a 
non-negative  integer  or  infinity. 

The  Cartesian  product  of  n  sets  X|,  .  .  .  ,Xn  is  denoted  by  either  Xj  X  XoX  •  ••  xX„  or 
X;,  and  defined  as 

n-X,  =X,  X  X,x  ...  X.x.  s  ((X, . X.)  ...x,  €  X,  ). 

1-1 

n 

If  x=(x,,  .  .  .  ,x„)€  nXj,  then  Xj^Xj  is  called  the  i-th  coordinate  of  x. 
i-i 

n 

Given  sets  Xi,...,Xn,  the  function  pj  nxi  -Xj,  where  1  ^  j  <  n,  defined  by 

i-l 

Pj(xi,  •  ■  •  ,Xj,  •  •  ■  ,Xn)  =  Xj,  is  called  the  projection  onto  the  j-th  coordinate. 

n 

Conversely,  given  a  function  q:X  — ►  11^ >  written  as  an  ordered  n- 

i-l 

tuple  of  functions  q  =(qi,  .  •  .  ,qn),  where  each  qj ;  X  — ►  X;  is  defined  by  q;  =spj  q;  that 
is,  qi(x)=Pi(q(x)).  The  function  qj  is  called  the  ith-coordinate  function  of  q. 

We  also  found  it  useful  to  employ  the  extended  real  number  system  in  both,  image  pro¬ 
cessing  theory  and  practice.  The  extended  real  number  system  Roo  consists  of  the  real 
number  system  R  to  which  two  symbols,  -HX3  and  — oo.  have  been  adjoined.  The  arithmetic 
and  logic  operations  of  R  are  extended  to  R^  as  follows:  Let  r€R,  then 
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(a) 

-oo  <  r  <  oo, 

r  V  oo  =  oo. 

r  A-oo  =  —00, 

(b) 

r  -H  oo  =  -hoo 

r  —  oo  =  — oo 

o 

II 

II 

"1 

+00  — oo 

(c) 

If  r  >  0,  then 

r  •  (+oo)  =  -HX) 

r  •  (— oo)  =  — oo 

(d) 

If  r  <  0,  then 

r  •  (-l-oo)  =  — oo 

r  •  (— oo)  =  -Hx) 

The  notation  V  and  A  in  (a)  above  denote  maximum  and  minimum,  respectively. 


2.  ALGEBRAS  AND  VALUE  SETS 

Intuitively,  an  algebra  is  simply  a  set  together  with  a  finite  number  of  operations  (rules) 
for  transforming  one  or  more  elements  of  the  set  into  another  element  of  the  set.  In  precise 
mathematical  terms,  an  algebra  is  a  system  {  F,0  }  in  which  F  is  a  non-empty  set  and  0  = 
{  o„  }  is  a  set  of  finitary  (n-ary)  operations,  where  each  operation  Oi^  G  0  is  a  mapping 

k 

0|, .  n  F  — ♦  F 

i-l 


For  example,  if  k=l  or  k=2,  then  is  a  unary  or  binary  operation,  respectively. 

Whenever  the  set  of  operations  0  is  tacitly  understood  it  is  customary  to  let  F  denote 
the  algebra  {  F,  0  }. 

Example  2.1. 

The  real  numbers  R  together  with  the  operations  of  addition,  multiplication,  and 
maximum  forms  an  algebra  that  falls  into  a  general  class  of  algebras  known  as  vec¬ 
tor  space  lattices. 

The  set  Zg  =  {  0,1, 2,.. ..,7  }  with  the  operation  of  addition  modulo  8  is  an  algebra. 
This  specific  algebra  belongs  to  a  class  of  algebras  known  as  cyclic  abelian  groups. 

n 

Suppose  F  =  nFj,  where  each  Fj  is  an  algebra  whose  set  of  finitar.y  operations  is  Oj. 

1—1 

The  set  of  naturally  induced  operations,  0,  on  F  is  defined  as: 

n 

0={  o€  n  Oj ;  the  coordinates  of  o  have  the  same  arity  }. 

1—1 

Thus,  if  o=(0|, . ,oJ  and  o-,  is  a  binary  operation  on  Fj  for  some  i  between  1  and  n,  then 


11 


o€  0  only  if  each  of  the  remaining  coordinates  oj  is  a  binary  operation. 

If  f =(fj, . ,fn)€F  g={gii . .gn)^F,  then  for  o=(oj, . ,o„)€  0  we  define 

fog=(f,o,g„ . ,f„o„g„) 

II 

The  set  F  =  n  F;  together  with  its  naturally  induced  set  of  operations  is  called  a  value  set. 
i-l 

It  follows  that  every  value  set  is  an  algebra. 

a 

If  F  =  n  F;  and  Fi=F:  for  all  i,j  =  l,2,....,n,  then  F  is  called  a  homogeneous  value  set, 
i«l  ^ 

Otherwise  F  is  called  a  heterogeneous  value  set.  If  F  is  homogeneous  and  o=(oi, . ,On)  is 

aji  operation  on  F  with  0;=Oj  for  all  i,j  =  l,2,....,n,  then  o  is  called  a  homogeneous  opera¬ 
tion,  otherwise  o  is  called  a  heterogeneous  operation.  Whenever  o  is  homogeneous,  then  it 
is  customary  to  use  o-,  to  denote  the  operation  o. 

Ebcample  2.2. 

For  i=l,2,  and  3,  let  Fj  =R  and  0|  =*{+,*, V}.  Then  F  =R^  and,  for 
example,  +=(+,+, +),  V =(V,\/,\/),  and  o=(+,*,V)  are  all  elements  of  0. 

Hence,  if  a3:(ai, 3.^,3^)  and  b=(bi,b2,b3)  are  elements  of  R^,  then 
a4-b=(aj+bi,a2-fb2,a3-l-b3),  aVb=«{aj Vbj,a2Vb2,a3Vb3),  and 
a  o  b  =(ai -|-bi,a2  *  b2,a3  Vbj). 

Some  of  the  more  common  examples  of  value  sets  used  in  image  processing  tasks  are  the 
sets  Z,  R,  C,  and  Zju  together  with  their  usual  arithmetic  and  lattice  operations. 

Given  a  value  set  {  F,0  },  then  in  addition  to  the  operations  given  by  0,  the  image 
algebra  of  F  valued  image  as  defined  in  subsequent  sections  also  incorporates  the  operations 
of  U,  n,  \  ,  choice  function,  and  cardinality  function  on  subsets  of  F,  where  \  denotes  set 
subtraction.  Thus,  if  F  =  R,  then  the  operations  are  the  usual  arithmetic  operations  +  ,  -  , 

*  ,  V,  A,  etc.,  together  with  union,  intersection,  subtraction,  choice  function,  and  cardinality 
function  on  subsets  of  R. 

3.  COORDINATE  SETS 

For  i  =  1,2, . ,k,  let  R;”'  denote  nj-dimensional  Euclidean  space,  each  having  a  possibly 

different  coordinate  system  assignment.  Thus,  for  example,  if  but  nj=nj=2,  then  R;" 
need  not  have  the  same  coordinate  system  as  Rj",  even  though  topologically  both  are 
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equivalent  to  R^.  Let  X^CRj”'  be  a  non-empty  subset  of  nj-dimensional  Euclidean  space.  A 
coordinate  set  X  is  a  set  of  form 

x=  nx;. 

i-i 

In  particular,  if  =:  a  and  k  =  1,  then  XCR". 

The  elements  of  X  will  be  denoted  by  bold  lower  case  letters.  Thus,  if  x€XCR",  then 
X  is  of  form  x=(xi,X2, .....  ,Xa),  where  each  coordinate  Xj  is  a  real  number. 

Ebcample  2.3. 

If  Xi  =  {(xi,X2) :  1  <Xi  <3,  1  <X2<4,  Xj,  Xo 6  Z  }  and  k=l,  then  X  =  Xj  C  R^  is 
the  array  shown  in  Figure  1. 
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Figure  1.  A  Finite  Rectangular  Array  in  Euclidean  2-Space  R^ 

If  all  the  factors  of  X  have  identical  underlying  coordinate  systems,  then  X  is  called  a 
homogeneous  coordinate  set,  otherwise  X  is  called  a  heterogeneous  coordinate  set.  In  Exam¬ 
ple  2.3  above,  X  is  a  homogeneous  coordinate  set.  We  now  provide  an  example  of  a  hetero¬ 
geneous  coordinate  set. 

Example  2.4. 

Let  X=XiX  Xj,  where  X^  and  Xo  denote  two  finite  rectangular  arrays  in 


13 


Z^—Zx  Z  having  different  underlying  coordinate  systems.  We  may  visualize  Xj 
and  X2  in  terms  of  cells  as  shown  in  Figure  2,  where  the  center  of  each  cell 
corresponds  to  a  point  in  the  respective  array  and  neighboring  cells  correspond  to 
neighboring  points.  Note  the  size  of  cells  in  one  array  differ  from  those  in  the 
other.  This  situation  could  correspond  to  two  sensors  of  different  resolving  power 
viewing  approximately  the  same  scene. 


Figure  2.  Two  Overlapping  Coordinate  Sets  of  Different  Resolution 

It  follows  from  the  definition  that  coordinate  sets  can  be  rectangular,  hexagonal,  toroidal 
discrete  arrays  as  well  as  stacks  of  rectangular  arrays  of  different  resolution  or  Euclidean 
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space  R°.  Providing  coordinate  sets  with  such  wide  varieties  of  shapes,  sizes  and  dimensions 
allows  to  a  coherent  approach  for  the  manipulation  of  differently  digitized  images  and 
multisensor/multidata  images. 

Image  algebra  operations  acting  on  coordinate  sets  are  operations  between  subsets  of 
coordinate  sets  as  well  as  operations  between  coordinate  points.  In  particular,  operations 
between  subsets  of  coordinate  sets  are  U,  D,  \ ,  choice  function,  and  cardinality  function. 
Operations  on  coordinate  sets  are  the  usual  operations  between  coordinate  points,  i.e.,  vector 
addition,  scalar  and  vector  multiplication,  dot  product,  etc. 

4.  IMAGES 

Thus  far  we  have  defined  two  types  of  objects,  value  sets  and  coordinate  sets.  These 
sets  and  their  elements  constitute  some  important  operands  of  the  image  algebra.  However, 
the  most  fundamental  of  the  algebra’s  operands  are  images.  The  most  general,  yet  useful, 
mathematical  definition  of  an  image  involves  the  previously  defined  concepts  of  value  sets 
and  coordinate  sets. 

Given  a  coordinate  and  value  sets  X  and  F,  respectively,  then  an  F  valued  image  a  on 
X  is  the  graph  of  a  function  a:X  — »  F.  Thus,  am  F  valued  image  a  on  X  is  of  the  form: 

a  =  {  (x,a(x)) :  x  €X  *(x)  €  F  } 

The  set  X  is  called  the  set  of  image  coordinates  of  a,  and  the  range  of  the  function  a 
(which  is  a  subset  of  F)  is  the  set  of  image  vadues  of  a.  If  the  vadue  set  F  is  the  set  of  posi¬ 
tive  real  numbers,  then  the  image  values  of  a  are  also  cadled  gray  values.  The  set  of  all  F 
valued  images  on  X  is  denoted  by  F^. 

If  the  value  set  F  =  R  or  F  =  R“,  then  we  are  dealing  with  real  valued  or  n- 
dimensional  vector  valued  images,  respectively.  Similarly,  replacing  F  by  Z,  C,  or  Z.*,  pro¬ 
vides  for  integral,  complex  or  finite  digital  images,  respectively.  Replacing  F  by  Rg^  proves 
useful  when  manipulating  raw  radar  images,  which  often  contain  spatial  locations  with  out- 
of-range  values  or  with  no  signal  values  (called  missing  values).  These  locations  can  be 
assigned  special  symbols  corresponding  to  oo  and  — oo,  respectively. 

The  examples  of  the  preceding  paragraph  should  make  it  clear  that  the  various  choices 
for  F  allow  for  a  far  greater  variety  of  image  operands  than  are  currently  used  by  the  image 
processing  community. 
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5.  BINARY  AND  UNARY  OPERATIONS  ON  IMAGES 

The  elementary  operations  between  F  valued  images  are  the  natural  induced  operations 
of  the  algebraic  system  F.  For  example,  for  real  valued  images  (i.e.  elements  of  R^),  the 
operations  are  the  elementary  operations  induced  by  the  vector  lattice  (=  a  vector  space 
which  is  also  a  lattice)  R.  Thus,  the  basic  real  valued  image  operations  reflect  the  arithmetic 
and  logic  operations  on  R.  In  particular,  the  binary  operations  of  addition,  multiplication, 
exponentiation  and  maximum  on  are  defined  as  follows: 

Let  a,  b  €  R^.  Then 

a  +  b  =  {  (x,c(x))  ;  c(x)  =  a(x)  +  b(x),  x  €  X}  (1) 

a  *  b  =  {  (x,c(x))  ;  c(x)  =  a(x)  *  b(x),  x  €  X}  (2) 

a  V  b  =  {  (x,c(x))  :  c{x)  =  a(x)  V  b(x),  x  6  X}  (3) 

These  are  the  basic  operations  for  real  valued  images.  As  complex  numbers  are  not 

endowed  with  a  natural  lattice  structure,  only  operations  1  and  2  are  basic  operations 
between  complex  valued  images.  With  the  exception  of  operation  6  below,  the  remaining 
binary  operations  on  real  or  complex  valued  images  could  be  obtained  from  these  basic 
operations  either  directly  or  in  terms  of  series  expansion. 

a**  s  {(x,c(x))  ;  c(x)  =»  a(x)‘^*)  if  a(x)  #  0,  else  c(x)  =  0,  x  €  X  } .  (4) 

We  restrict  this  binary  operation  to  those  pairs  of  images  a,b  for  which  a(x)*^’'^€R 
whenever  a(x)7^0.  The  inverse  of  exponentiation  is  defined  in  the  usual  way  by  taking  the 
logarithm.  In  particular,  we  define: 

logab  =  {(x.c(x)) :  c(x)  =  log^,j)b(x),  x  G  X  } .  (.5) 

As  for  real  numbers,  log^b  is  defined  only  for  those  images  a  and  b  for  which  a(x)  >  0 
and  b(x)  >  0  for  adl  x  €  X  The  next  basic  binary  operation,  called  the  dot  product,  distin¬ 
guishes  itself  from  the  above  five  in  that  its  output  is  not  an  image  but  a  real  number.  Let  X 
be  finite,  then  the  dot  product  is  defined  as: 

a«b  =  S  a(x)b(x)  (6) 
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An  image  a  is  called  a  constant  image  if  the  same  value  is  assigned  by  a  to  each  element 
of  its  coordinate  set;  i.e.  if  a(x)  =  k  for  some  real  number  k  and  for  all  x  €  X. 

Two  constant  images  that  are  of  utmost  importance  in  the  Image  Algebra  are  the  zero 
image,  defined  by  0  =  {(x,0):  x€X  },  and  the  unit  image,  defined  by  I  — {(x.l);x€X}. 

Suppose  k  G  R  and  a  is  a  constant  image  with  a(x)  =  k.  Then  we  define: 
b''  =  b*  and  k**  =  a'’, 
kb  =  a*b  and  k+b  =  a  +  b, 
logjjb  =  log^b,  of  course  k>0  and  b(x)>0  for  all  x. 

We  note  that  exponentiation  is  defined  even  when  a(x)  =0.  Subtraction,  division  and 
minimum  are  defined  in  terms  of  the  basic  operations  and  inverses.  Specifically: 

a  -  b  =  a  +  (— b),  where  -b  =  {(x,-b(x));  (x,b(x))  G  b  }, 
a/b  =  a*b~\  and 
a  A  b  =  —(—a  V  -b). 

The  images  0  and  1  have  the  obvious  property  a+O  =  a  and  a*l  =  a.  On  the  other 
hand,  b*b~^  does  not  necessarily  equal  1.  However,  b*b~^*b  =  b.  For  this  reason  b~^  is 
called  the  pseudo  inverse  of  b.  Inequalities  between  images  are  defined  in  terms  of  maximum 
and  minimum.  Thus,  for  example,  a<b  if  aVb=b.  These  observations  show  that  the  ring 
(R^,+,*)  and  the  lattice  (R^,V.<)  behave  very  much  like  the  ring  and  lattice  of  real 
numbers. 

There  are  various  useful  unary  operations  definable  in  terms  of  the  basic  binary  opera¬ 
tions.  For  example,  we  have  already  provided  the  definitions  of  k**  and  logj^b,  the  exponen¬ 
tial  of  an  image  and  the  logarithm  of  an  image  b  to  the  base  k,  respectively.  In  particular, 
the  exponential  of  an  image  b,  exp(b)  =6**  =  a'’,  and  the  natural  logarithm  of  b  is  defined 
as  lnb=log,b,  respectively,  where  a  is  the  constant  image  defined  by  a(x)  =  e  for  all  xGX 

If  a  and  b  are  images,  then  the  characteristic  value  Y>b(*)  ^  defined  as 
X>b(a)  =  [{a-b)  V  Oj"'  *  l(a-b)  V  0| 

Thus, 

X>b{a)  =  {  fx.c(x))  :  c(x)=l  if  a(x)>b(x),  else  c(x)=0  }. 
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Obviously,  if  a=b,  then  (a  -  b)  V  0  =  0  and  X>b  (*^)  “  0“**0  =  0  since  by  definition  of 
exponentiation,  0~^  =  0.  The  complement  of  an  image  a  is  denoted  by  a  and  is  defined  as 
a  =  1  —  a*a~^  The  remaining  characteristic  functions  of  images  can  be  defined  in  a  similar 
fashion,  using  complementation  and  products.  For  example, 

=  X>b(a) 
and 

>(b(a)  =  X^(a)  * 

Whenever  b  is  the  constant  image  with  gray  values  equal  to  k  it  is  customary  to  replace  b 
by  k  in  the  above  definitions.  Figure  3  below  provides  an  example  of  the  operation  X>it(a), 
where  k  =  15. 


Figure  3.  Example  of  the  Operation  ,\>is(a) 

Observe  that  the  characteristic  function  provides  a  good  example  as  to  the  use  of  the 
inverse  of  an  image  and  the  use  of  pointwise  multiplication  and  maximum  operations. 
.\nother  similar  example  is  the  absolute  value  of  an  image  which  is  defined  as  1  al  =aV(— a). 
These  examples  show  that  there  are  various  image  operations  of  different  degrees  of  complex¬ 
ity  that  can  be  defined  in  terms  of  the  more  basic  first  six  operations  listed  above.  Thus, 
instead  of  listing  finite  strings  of  basic  operations  in  order  to  represent  the  characteristic 
function  of  an  image  a  Avith  respect  to  any  non-empty  set  SCR,  Ave  use  the  more  economical 
symbols,  ?<s(a)  and  lal  ,  in  order  to  represent  the  image  { (x,c(x))  :  c(x)=l  if  a(x)€S.  else 
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c(x)  =  0  }  and  aV(-a),  respectively. 

Thus  far  unary  operations  have  been  defined  in  terms  of  the  basic  binary  operations  and 
we  observed  that  they  aire  not  elemental  in  the  sense  of  defining  the  Image  Algebra.  The 
next  set  of  unary  operations  could  also  be  defined  by  the  above  named  basic  operations  in 
terms  of  series  approximation.  However,  this  would  defeat  the  goal  of  providing  a  simple 
lamguage  for  image  processing  tasks. 

The  basic  unary  operations  on  are  the  functions  available  in  most  high  level  pro¬ 
gramming  languages.  In  fact,  any  function  /  ;  R  — ►  R  induces  a  function  R^  — ►  R^,  again 
denoted  by  /,  and  defined  by 

/(a)  =  {(x.cW):c(x)  =/{a(x))} 

For  example,  sin(a)  =  {(x,sin(a(x))):  x  €  X  }. 

All  operations  defined  in  this  section  with  the  exception  of  those  involving  the  lattice 
operations  V  and  A  and  relations  <  and  ^  can  also  be  applied  to  complex  valued  images, 
i.e.  elements  of  C^.  Similarly,  as  mentioned  before,  the  operations  on  are  the  natural 
induced  operations  of  the  algebraic  system  F.  For  instance,  if  F  =R“,  then  a(x)GR"  is  of 

form  ai(x)={ai(x),a2(x) . ,an(x)),  with  a^(x)€R  for  is=l,2,....,n,  and  all  previously  defined 

operations  are  defined  componentwise.  In  particular, 

a(x)  -l-b(x)  =(ai(x)  +bj(x) . ,a„(x)  +b„{x)), 

a(x)  •  b(x)  =(ai(x)  •  bj(x), . ,a„(x)  •  b„{x)),  and 

sin(a(x))=<sin(ai(x)), . ,sin(a„(x))). 

When  considering  an  expression  of  the  form  sin(a),  one  rarely  thinks  of  the  image  a  as  a 
function.  However,  in  this  paper,  images  are  defined  as  functions,  namely  elements  of  F^, 
and  two  important  mathematical  notions  used  in  image  processing  are  the  restriction  and 
extension  of  a  function.  We  express  these  notions  as  two  basic  operations  of  the  image  alge¬ 
bra. 

Let  a  be  an  image  on  X  The  restriction  of  a  to  a  subset  Y  of  X  will  be  denoted  by 
a]  Y-  Here  a  user  would  specify  the  coordinate  set  Y  C  X  As  an  example,  he  could  set  Y  = 
(x€X  :  5  <1x1  ^0  }.  The  output  image  aiy  is  an  image  on  Y. 

Let  a  be  an  image  on  X  b  an  image  on  Y,  and  X  C  Y.  The  extension  of  a  to  b  on  Y 
is  defined  by 
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J(b.YV,>  if ’■ex 

ifxGY\X 


where  the  user  specifies  the  function  (image)  b  on  Y. 

Essential  to  the  definition  of  a  function  are  its  domain  and  range.  These  form  the 
remaining  basic  operations  on  images.  The  domain  function  of  an  image  a  is  defined  simply 
as 

Domain(a)  =  set  over  which  a  is  defined, 

e.g.,  Domain(al  y)  =  Y.  The  range  function  of  an  image  lies  on  the  other  side  of  the  spec¬ 
trum,  it  provides  the  set  of  values  assigned  to  coordinates  by  the  image.  Thus,  if  a  is  an 
image,  then 

Range(a)  =  set  of  values  determined  by  a, 
e.g.,  Range(al  y)  is  the  set  of  all  values  a  assumes  on  Y. 

6.  GENERALIZED  TEMPLATES 

In  terms  of  image  processing  applications,  templates  and  template  operations  are  the 
most  powerful  tool  of  the  Image  Algebra.  Generalized  templates  are  the  most  abstract 
objects  in  the  algebra.  One  reason  for  the  abstraction  is  that  our  definition  of  template 
unifies  and  generalizes  the  usual  concepts  of  templates,  masks,  windows  and  neighborhood 
functions  of  other  image  processing  formalisms  into  one  general  mathematical  entity. 

Before  defining  generalized  templates,  we  will  look  at  a  function  that  assigns  to  each 
point  of  a  set  a  pair  of  objects.  Henceforth,  let  X  and  Y  be  coordinate  sets,  F  a  value  set, 
and  let  2^  denote  the  power  set  of  X  is,  2^  is  the  set  of  all  subsets  of  X  Consider  a 
function  of  type 

t={ti,t2)  ;Y  -••2^xF^ 

Note  that  for  each  point  y  GY,  t(y)  =(ti(y),to(y))  where  ti(y)G2^  and  t2(y)  G  F^. 
That  is,  tj(y)  C  X  and  t2(y)  is  an  image  on  X 

For  subsequent  notational  convenience  we  will  denote  the  coordinates  of  this  type  of 
function  in  upper  case  script  and  lower  case  italics.  In  particular,  the  coordinates  of 
t=(ti,t2),  tj  and  to,  will  be  denoted  by  T  and  t,  respectively.  In  addition,  we  define 
ty  =  t(y).  Thus,  ty  is  the  image  ty  =  {(x,ty(x)):  xGX  }. 
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Henceforth,  all  value  sets  under  consideration  are  assumed  to  have  an  operation  of  addi¬ 
tion  and  an  additive  identity.  A  generalized  F  valued  template  t  from  Y  to  X  is  a  function 
t:Y— ►  2^xF^  whose  coordinates  (T,  t)  satisfy  the  property  that  for  each 
y€Y,  ^(x)=0  whenever  x^T^y).  Here,  0  denotes  the  additive  identity  of  F. 

An  equivalent  way  of  stating  this  definition  is  to  say  that  an  F  valued  template  t  from 
Y  to  X  is  a  pair  t  =  (T,  t),  where 

1.  T  Y  2*“;  i.e.,  T(y)  C  X  for  each  y €Y  and 

2.  t  Y  F*  with  ty  =  {(x,^(x)):  ty(x)=0  if  x$f  T(y),  x€X  }. 

The  set  T(y)  is  called  the  source  configuration  of  y;  y  is  called  the  target  point  of  T(y); 
and  the  values  <y(x)  for  x  €  T(y)  are  called  the  weights  of  t.  Figure  4  illustrates  these  con¬ 
cepts. 


Figure  4.  A  Pictorial  Example  of  a  Template  from  Y  to  X 

The  following  example  should  further  clarify  the  template  concept.  Let  XCR“  be  a 
rectangular  array,  let  y  =(x,y)  be  an  arbitrary  point  of  X,  Xj  =  (x,y— 1),  Xn  =  (x-(-l,y),  and 
X3  =  (x+l,y— 1).  We  now  define  a  template  t  =  (T,t)  from  X  to  X  by  defining  -  for  each 
y€X-  its  weights  as  ty(x)  =0,  ty(xi)  =  I,  ty(x2)  =2,  ^(xs)  =3,  and  ty(x)  =0  if  x  is  not  an 
element  of  the  source  configuration  T(y)  =  {y,Xi,Xo,X3}.  Thus  t  has  configuration  and 
weights  as  shown; 
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The  shaded  cell  in  the  pictorial  representation  of  t  indicates  the  location  of  the  target  point 

y 

Templates  are  used  to  define  those  image  transformations  which  make  use  of  all  image 
values  within  some  predescribed  configuration  of  the  source  or  input  image.  For  example, 
given  an  image  a  on  X  and  t  =(7,t)  a  template  from  Y  to  X  where  Y  may  be  of  an 
entirely  difierent  shape,  size  or  dimension  than  X  l^hen  an  operation  between  a  and  t  will 
transform  a  into  an  image  b  on  Y  where  each  new  pbcel  value  b(y)  is  computed  in  terms  of 
some  arithmetic  and/or  logic  combination  of  the  values  a(x)  and  <^(x),  where  x€  T(y).  Ini¬ 
tially,  however,  it  may  be  convenient  to  view  templates  as  masks  such  as  the  edge  masks 
used  in  the  Sobel  or  Kirsch  edge  detection  schemes.  The  template  configuration  correspond¬ 
ing  to  a  Sobel  or  Kirsch  edge  mask  is  always  a  3x3  neighborhood  configuration  about  the 
center  pixel.  However,  the  weight  functions  corresponding  to  the  two  Sobel  edge  masks  (or 
eight  Kirsch  edge  masks)  are  distinct.  Thus,  the  two  Sobel  masks  given  in  Example  2.6  are 
examples  of  two  distinct  templates  having  the  same  configuration.  It  is  important  to  realize, 
however,  that  the  notion  of  a  template  is  not  analogous  to  that  of  a  mask.  The 
magnification  and  Fourier  templates  presented  in  Example  2.10  and  Section  HI,  respectively, 
are  examples  of  templates  which  are  not  masks.  Here  the  weights  or  configuration  changes 
as  a  function  of  the  position  of  the  target  pixel.  In  addition,  the  target  pbcel  of  these  tem¬ 
plates  does  not  correspond  to  the  geometric  center  of  the  configuration.  In  general,  given 
two  distinct  templates  t  =  (T,t)  and  s  =  (S,s)  from  Y  to  X  one,  and  only  one,  of  the  follow¬ 
ing  three  cases  is  possible:  (i)  T and  t^s,  or  (ii)  7=5  and  or  (iii)  7 and  t=s. 

The  set  of  all  F  valued  templates  from  Y  to  X  will,  henceforth,  be  denoted  by  Fy_x- 

7.  OPERATIONS  BETWEEN  IMAGES  AND  TEMPLATES 

There  are  three  basic  template  operations  that  are  used  to  transform  a  real  valued 
image.  They  are  denoted  @,  <2l ,  and  Q  ,  and  called  generalized  convolution,  multiplicative 
maximum,  and  additive  maximum,  respectively.  For  complex  valued  images  only  one 
image-template  operation  is  defined,  namely  ©. 

A  template  operation  on  image  a  and  template  t  computes  a  pixel  value  c(y)  by  per¬ 
forming  the  basic  operation  of  addition  or  maximum  on  a  weighted  collection  of  all  pixel 
values  a(x)  with  coordinates  x  €  7(y).  In  particular,  if  a  G  and  t  G  Ry-^w  >  where 
X  and  W  are  subsets  of  the  same  Euclidean  space,  then 
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»®t  s  {(y,c(y)):  c(y)  =  E  a(x)-  ^(x),  y  €  Y}, 

xeT(y)nX 

=  {(y.c(y)):  c(y)  =  V  a<x)  <y(x),  y  €  Y},  and 

x€T(y)nX 

ants  {(y,c(y));  c(y)  =  V  a<x)+t^(x),  y  G  Y}, 

x€T(y)nX 

where  E  a(x)-^(x)  =0  whenever  T(y)nX  —  (t).  The  multiplicative  and  additive  max- 
x6Tty)nx 

imums  are  defined  only  if  T(y)nX  ^  (t>.  However,  if  we  extend  R  to  ,  then  ®  and  □ 

can  be  defined  by  V  a(x)-ty(x)  =  — oo  and  V  a(x)+^(x)  =  — oo,  respectively, 
x€T(y)nX  x€T(y)nX 

whenever  T^y)  flX  =  (2). 

The  complementary  operations  of  multiplicative  minimum  and  additive  minimum  are 
defined  in  terms  of  ®  and  Q  as  follows; 

a®t  s  — (a®  — t) 
a  Q  t  s  — { — a  n  — t) 

In  the  above  definitions  we  assume  that  71(y)nX  is  finite  for  each  y  GY.  However,  the 
definitions  extend  to  continuous  functions  a(x)  and  ty  on  compact  sets  T{y)  and  X  with  the 
exception  that  in  the  formulation  of  a0t  the  sum  is  replaced  by  an  integral.  That  is, 

c(y)  =  /  a(x)ty(x)dx 

T(y)nX 

Note  that  a  G  R^,  while  a@t  G  R^.  Thus,  template  operations  may  be  used  for  chang¬ 
ing  the  dimensionality  or  size  and  shape  of  images.  In  particular,  in  addition  to  the  usual 
local  or  global  convolutions — as  occur  in  edge  enhancement,  local  smoothing,  morphological 
operations  and  Fourier  like  transformations — template  operations  also  provide  a  tool  for 
image  rotation,  zooming,  image  reduction,  masked  extraction,  and  matrix  multiplication. 

The  operations  can  also  be  extended  to  lattice  ordered  vector  spaces.  If  Fj  is  a  vector  space 
over  the  field  Fg,  t  an  Fj  valued  template  from  Y  to  W  and  a  an  F2  valued  image  on  X 
then  a@t  is  an  F^  valued  image  (i.e.  a  vector  valued  image)  on  Y.  The  roles  of  a  and  t  can 
of  course  be  reversed,  that  is,  a  can  be  an  Fj  valued  image  and  t  an  Fo  valued  template. 
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We  note  that  other  logic/arithmetic  operations  between  image  values  and  template 
values  could  have  been  defined.  For  instance,  we  could  have  defined 

b(y)=  E  a(x)Vty(x). 
x€Tl(y)nX 

The  reason  for  not  defining  the  remaining  logic/arithmetic  operations  is  that  some  of  them 
are  easily  expressible  in  terms  of  the  three  basic  operations  defined  above,  and  no  application 
examples  could  be  found  for  the  remaining  combinations.  This  does  not  imply  that  they 
should  henceforth  be  ignored.  In  fact,  they  should  be  kept  in  mind  whenever  new  ima^e  pro¬ 
cessing  techniques  are  being  discovered  that  are  not  easily  expressible  within  the  current 
algebraic  structure. 

At  first  glance,  the  basic  image/template  operations  may  appear  somewhat  difficult. 
However,  the  examples  provided  below  and  those  given  in  References  10,5,11,  and  12  illus¬ 
trate  the  inherent  simplicity  and  power  of  these  operations. 

Example  2.5.  Local  Averaging 

Let  a  be  an  image  on  a  rectangular  array  X  C  R^,  Y  =  Z^,  and  t  G  Ry->y 
3x3  neighborhood  template  defined  as  follows; 


t(y)  = 


1 

1 

1 

1 

/  \ 

1 

s  / 

1 

1 

1 

1 

Then  — a@t  represents  the  image  obtained  from  a  by  local  averaging  since  the  new 
9 

pixel  value  is  given  by  c(y)=-^  S  a(x). 

9  x€r(y)nX 

Note  that  the  image  a© t  is  an  image  on  all  of  T?  with  zero  values  outside  of  the  array 
XC  T?.  Obviously,  computers  are  not  capable  of  storing  images  defined  on  infinite  arrays. 

Furthermore,  in  practice  one  is  only  interested  in  the  image  ■^(a©t)  restricted  to  the  array 
X  that  is  "(aStjlx,  where  1^  denotes  the  restriction  to  X. 
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This  problem  could  be  solved  as  follows:  Since  {a©t)lx  =  a®  (tlx);  <ihe  template 
8€Rx_*xi  defined  by  a  =  tlx  provides  for  the  desired  finite  image  -^aSs).  Thus,  the  ques¬ 
tion  arises:  "Why  not  simply  define  t  as  a  template  from  X  to  X  instead  from  T?  to  Z^?" 

The  reason  for  defining  the  template  as  we  did  is  that  this  template  can  be  used  for 
smoothing  any  2-dimensional  image  independent  of  its  array  size  X,  because  when  defining 
an  image  b  in  a  program  one  must  declare  its  dimensions,  i.e.  the  size  of  its  underlying  array 
(number  of  addresses)  Y.  Hence,  the  program  statement  b  =  a®  t  is  equivalent  to  the  image 
aJge bra  statement  b  =a@(tlY),  where  in  a  program  the  equal  sign  "=”  means  "replace"  b  by 
a©t.  Thus,  a  programmer  is  not  faced  with  the  task  of  redefining  t  for  a  different  sized 
image,  as  would  be  the  case  if  he  had  defined  t€Rx_»x  ^  given  X 

Example  2.6.  Sobel  Edge  Detection 

Let  a  be  an  image  on  a  rectangular  array  X  C  T?,  and  s^,  s^,  the  templates  shown 
below  defined  on  all  of  Z^.  The  image  algebra  expression 

l(a©8j‘^-l-(a®8y)Y^“, 

where 


-1 

0 

1 

.2 

/■  V 
0 

S  / 

2 

-1 

0 

1 

-1 

.2 

-1 

/  V 

0 

0 

0 

^  y 

1 

2 

1 

represents  the  Sobel  edge  enhanced  image.  Figure  5  provides  an  example  of  edge 
detection  using  the  Sobel  algorithm. 
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Figure  5.  input  image  (left)  and  Sobel  Edge  Enhanced  Image  (right) 


Example  2.7.  Geometric  Edge  Detection 

Geometric  edge  transforms  are  edge  enhancement  transforms  that  empioy  the 
image  transforms  a®  t  or  a  Q  t.  They  do  not — as  compared  to  the  Sobel 
Transform — make  use  of  the  notion  of  "derivative."  The  following  example  of  a 
geometric  edge  transform  mimics  the  Sobel.  Its  algebraic  formulation  is  as  follows; 

[(a  ®  tj  —  a  ®  to)^  +  (a  ®  tg  —  a  ® 

where 


1 

2 

0 

1 

1 

7 

0 

2 

_ 

1 
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Figure  8.  Example  of  a  Geometric  Edge  Transformation 
Example  2.8.  Dilations  and  Erosions 

We  present  this  example  for  readers  familiar  with  the  two  basic  notions  of  dilation 
and  erosion  that  define  all  image  processing  schemes  based  on  mathematical  mor¬ 
phology.  A  three  dimensional  dilation  of  an  image  a  by  a  structuring  element  t 
can  be  expressed  as  a  Q  t.  As  before,  the  corresponding  minimum  operation  Q 
can  be  defined  in  terms  of  0  by  a  Q  t  =  -(-a  0  -t).  In  particular,  if  t  has  only 
zero  weights  and  the  standard  3x3  neighborhood  configuration,  then  a  Q  t 
denotes  the  image  obtained  from  a  by  replacing  each  pixel  of  a  by  the  minimum  of 
the  values  in  its  immediate  neighborhood.  Figure  7(b)  represents  the  dilation  a0  t 
of  the  input  image  a  shown  in  Figure  7(a),  while  Figure  7(c)  represents  the  erosion 
of  the  image  a  0  t,  namely,  (a  0 1)  Q  t.  The  template  t  used  in  this  example  is 
defined  as  follows: 


I 
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Figure  7.  Template  t  of  Example  2.8 


(a)  The  Input  Image  a  (b)  The  Dilated  Image  b==a  Q  t  (c)  The  Eroded  Image  b  Q  t 


Figure  8.  Example  of  Dilation  and  Erosion 
Three  dimensional  erosion  of  an  image  a  by  a  structuring  element  t  is  given  by  the 
expression  a  Q  -t,  since  each  new  pixel  value  is  given  by 

c(x)  =  A  (a(y)  -  t,(y))  =  min{a^y)  -  t,(y) ;  y  E  T(y) }. 

y6T(x) 

It  is  important  to  note  that  a  template  may  vary  at  different  locations  in  both  shape  and 
weights.  Thus  the  expression  a  (3  t  may  represent  a  far  more  complex  algorithm  than  a 
simple  dilation.  In  short,  mathematical  morphology,  as  used  in  actual  image  processing,  is  a 
special  minor  substructure  of  the  mathematical  structure  represented  by  the  Image  Algebra, 
Reference  13. 


Before  discussing  parameterized  templates,  we  introduce  two  special  but  very  important 
operations  which  are  defined  in  terms  of  the  elementary  operations. 

The  sum  of  an  image  a  on  X  is  defined  as 

Sa  =la*l 

and  the  maximum  of  a  as 

Va  =  £(a©t), 

where  t  is  a  template  from  Y  =  {0}  to  X  defined  by  7(0)  =  X  and  fy(x)  =  1  for  each  x  €  X 
Note  that  Sa  =  S  a(x)  is  a  real  number  and  a©t  is  an  image  consisting  of  a  single  point. 

An  example  that  exhibits  the  brevity  of  image  algebra  code  and  involves  image  summa¬ 
tion  is  the  Euler  number  of  a  Boolean  image. 

Example  2.9.  Euler  Number 

The  Euler  number,  E,  of  a  Boolean  image  is  the  difference  between  the  number  of 
objects  (connected  components  having  pixel  value  1)  and  the  number  of  holes  in  the 
objects.  It  can  be  computed  locally,  thus  in  real  time  on  massively  parallel  archi¬ 
tectures.  The  algorithm  for  computing  E  is  especially  short  in  terms  of  the  Image 
Algebra: 

E  =  E[x2(a  ©  t)  -  X7(a  ©  t)  -+•  Xio(a  ©  <^)|. 
where  a  denotes  the  input  image  and  t  the  template 

t(y)  = 

8.  PARAMETERIZED  TEMPLATES 

Let  X,  Y  be  coordinate  sets,  F  a  value  set  with  identity,  and  P  a  non-empty  set.  .\ 
parameterized  F  valued  template  from  Y  to  X  with  parameters  in  P  is  a  function  of  form 

t;  Yx  P  -*•  F^ 


with  t(y,p)  =(Tl(y,p),^p)  having  the  property  that  lor  each  (y,p)€Yx  P  ,^p(x)  =0  when¬ 
ever  x^71[y,p).  Here  again  we  use  the  notation  ^  p  =  <(y,p). 

The  set  P  is  called  the  set  of  parsuneters  and  each  p€P  is  called  a  parameter  for  t. 

Given  a  parameterized  template  t:  Yx  P  — ►  2^X  F^,  we  can  define  a  family  of  F 
valued  templates  from  Y  to  X: 

tpiY  — *'2^X  F^ 

by  defining  for  each  p  G  P 

*p(y)  =  %.p)- 

Thus,  a  parameterized  F  valued  template  from  Y  to  X  gives  rise  to  a  family  of  regular 
F  valued  templates  from  Y  to  X  namely  { tp  :  pGP  }.  The  following  two  examples  should 
help  clarify  these  notions. 

Example  2.10.  Image  Magnification 

Suppose  X  C  R~  is  an  m  X  n  array,  Y  =  Z",  P  =  {  p  ;  p={xo,k),  where  x^  G  X.  k  a 
positive  integer  },  and  a  an  image  on  X.  Given  a  pair  of  real  numbers  r  =  (ri,ro), 
define  [r]  =  ([ri],[r2]),  where  [rj]  denotes  truncation  of  rj  to  the  nearest  integer.  For 
each  y  G  Y  and  p=Kxo,k),  let  Tp(y)  =  {  x  :  x  =  [(y  -  Xo)/k  +Xo]  }  and  <^_p(x)  =  1  if 
X  G  Tp(y).  Then  b  =  a  ©  tp  represents  the  magnification  of  a  by  the  factor  k 
about  the  point  Xq.  Thus,  once  this  parameterized  template  has  been  defined,  all  a 
potential  user  of  this  template  needs  to  supply  is  the  magnification  factor  k,  the 
point  about  which  to  magnify  the  image,  and  -  in  order  to  retain  all  the  information 
-  declare  b  to  be  of  at  least  dimension  kmX  kn.  This  example  also  shows  how  a 
template  transformation  is  capable  of  changing  the  size  of  an  image.  In  the  figure 
below,  the  image  on  the  right  represents  the  magnification  of  the  image  on  the  left 
by  a  factor  of  2. 
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Figure  9.  An  Example  of  Image  Magnification  by  a  Factor  of  2 
Exampl'i  2.11.  Image  Rotation 

Suppose  XCR^  is  an  mX  n  array,  Y  =  Z^,  P  =  {p:p=(xo,^),  where  x^GX,  9  an 
angle  between  0  and  27r  },  and  a  an  image  on  X.  Declare  c  to  be  of  sufficiently 
large  dimension  such  that  when  a  is  rotated  through  an  angle  9  about  a  point  Xq, 
the  rotated  image  fits  into  the  array  determined  by  c.  We  view  the  points  of  X 
and  Y  as  cells,  i.e.,  a  unit  box  with  center  at  each  point.  For  each  y€Y  and  p  = 
(xq,^),  let  the  configuration  Tp(y)  be  the  set  of  cells  of  X  having  non-empty  inter¬ 
section  with  the  cell  z,  where  z  is  obtained  from  y  after  rotating  y  clockwise 
through  the  angle  9  about  the  point  Xg,  as  shown  in  Figure  9. 
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Figure  10.  Target  Point  y  and  Source  Configuration  7(y)  of  the  Rotation  Template  t 

The  weights  ty, p(xi),  where  Xj  €  Tp(y),  are  defined  by  ly,p(xi)  =  area  of  intersection 
of  the  preimage  of  y  with  Xj.  The  rotated  image  c  is  now  given  by  c  =  a©t,  where 

c(y)  =  ,a(x)‘y,p(x)- 

*€Tp(y) 

Figure  11  provides  an  example  of  image  rotation  using  the  template  just  defined. 
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Figure  11.  Rotating  an  Image  Through  15,  30,  45,  60,  and  75  Degrees 


Example  2.12.  The  Kirsch  Edge  Detector 

The  standard  formulation  of  the  Kirsch  edge  detection  algorithm  is  to  replace  each 
pixel  a(y)  of  the  input  image  a  by 

c(y)  =  max{l,  max{  I  5(24  +  aj+j  +  ai^.^)  -  3(0.;^^  +  ...  +  3^+7)  1 ;  i  =  1,2,  ...,  8  }}, 

where  the  addition  of  subscripts  is  mod  8  and  the  a^’s  denote  the  following  eight 
neighbors  of  a(y): 


^3 

^5 

/"  'X 
a(y) 

ai 

^6 
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Figure  12.  Example  of  the  Kirsch  Edge  Detection  Algorithm 

9.  OPERATIONS  BETWEEN  TEMPLATES 

The  operations  ©,  ® ,  and  0  between  images  and  templat.iS  generalize  to  operations 
between  templates.  In  particular,  if  t  is  a  real  valued  template  from  Y  to  X  and  s  is  a  real 
valued  template  from  X  to  W,  then  we  define  a  new  template  r  =8  ©  t  from  Y  to  W  by 
defining  its  weight  function  r  by 
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ifwe^(y) 

•  x67ty)  ' 

“  0  otherwise 

where  w€W  and  ^(y)  *=  U  S(x). 

*€T(y) 

Similarly,  for  r  =8  Q  t  we  define  r  by 

V{  ty(x} + 8,(w) :  X  €  T(y )  and  w  €  S{x)} 
ry(w)=  ifw€-e(y) 

0  otherwise 

amd  for  r  =8  ®  t  we  define  r  by 

V{  iy(x)«^(w) :  X  €  T{y)  and  w  G  5(x)} 
ry(’w)=  ifw€^(y) 

0  otherwise 

The  complementary  operations  E3  and  ®  are  defined  by 

8  £3  fc  =  — ( — 8  Q  — t) 

8  ®  t  =  —(8®  -t) 

where  wGW  and  in  each  case  R,{y)  =  U  5(x). 

xeT(y) 

If  8  and  t  are  templates  from  Y  to  X,  then  addition,  multiplication,  and  maximum 
between  s  and  t  are  defined  pointwise.  In  particular,  we  define 

8  +t  by  {«  +  t)y  =«y+ty 

8  *  t  by  («  *  t)y  =  8y  *  ty 

and  8  V  t  by  («  V  t)y  =  fly  V  ty 

The  source  configuration  of  each  of  these  templates  at  the  target  point  y  is  given  by 

S{y)  U  T(y). 

Initially,  these  definitions  seem  to  be  fairly  complex.  The  following  examples  serve  to 
clarify  these  definitions  and  should  provide  a  better  understanding  as  to  how  composition  of 
templates  is  accomplished. 
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Ebcample  2.13. 


1 

/  V 

Let  s  = 

1 

2 

1 

t  = 

3 

-1 

then 


1 

1 

/  \ 
5 

\  / 

1 

-1 

0 

0 

/  \ 
6 

\  y 

0 

0 

1 

1 

/  \ 
3 

1 

0 

1 

2 

1 

3 

6 

3 

-1 

-2 

-1 

As  can  be  seen,  template  addition  and  multiplication  do  not  pose  any  great  conceptual 
difficulties.  In  contrast,  the  operations  s  @  t,  s  ®  t,  and  s  0  t  are  somewhat  more  intri¬ 
cate.  For  this  reason  we  provide  a  step  by  step  analysis  of  these  operations.  We  start  by 

computing  the  configuration  ^(y)  =  U  S(x)  at  a  point  y€Y.  Specifically,  if  y€Y,  then 

x€r(y) 
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T(y)  is  of  the  form 


7’(y)  = 


where  Xj,  x^,  Xg  €  X  with  y  =  x<,.  Hence  S(x^),  S(x^),  and  <?(xg)  are  located  as  follows. 


^(xi) 

^(Xo) 


7 - ^ 

/  \ 
Xo 

S  “/ 

/  \ 

and  the  configuration  is  the  union: 


^(y)  *=  ^ gUyj  S{x)  *  Six,)  U  S(x.,)  U  Six,) 


We  now  compute 


r  (w)  = 


>M*)*x(w), 

x6T(y) 

0  otherwise 


if  weje(y) 


Since  ry(w)=0  if  (w)^.S(y),  we  need  only  compute  ry(w)  for  w€.C(y).  In  order  to  illustrate 
this  task,  we  relabel  the  points  of  JZiy)  as  follows: 


Wg 

< 

CO 

7  \ 

^  V 

^6 

^7 

< 

00 

^9 

Here  Wo  =  Xj,  Wo  =  Xj,  W5  =  y  =  Xn  and  Wg  =  Xg.  Computation  of  ry(wi)  for  each  i  =  1,  .... 
9  now  follows  from  the  definition: 
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A(*K(Wi),  where  W;  6  S{x) 
x€T(y)  ' 

Specifically, 

=  <y(*iK,(wi)  =  1*1  =  1 
Tylwj)  =  ty(xi)«,^(w2)  =  1*2  =  2 
^yM  =  *y(XlK,(W3)  =  1*1  =  1 

»-y(w4)  =  =  3*1  =  3 

’•yC'Ws)  =  *y(x2)« Jws)  =  3*2  =  6 
»‘y(w8)  =  ^W^xi'We)  =  3*1  =  3 

ry(w7)  =  Iy(x3)«^(w7)  =  -1*1  =  -1 

»‘y('»»'8)  =  ty{^z)»xS'"&)  =  “1*2  =  -2 

^(^9)  =  <y(x3)*x,(W9)  =  -1*1  =  -1 

Note  that  since  Wi^S(x2)  U  S(x3),  the  definition  does  not  require  that  the  quantities  ex„(wj) 
and  enter  into  the  computation  of  Wj.  However,  since  and  nre  both 

zero,  there  is  no  harm  in  including  them  in  the  sum.  In  fact  for  the  operation  @  it  is  always 
possible  to  compute  Wj  by  relaxing  the  condition  stated  below  the  S  and  including  extra 
terms,  all  of  which  equal  zero.  Thus 

3 

=  2<y(Xj)«x,(wi) 

=  <y(*l)«x,('»^i)  +  *y(x2)«x„{Wi)  +  ly(X3)«x,(Wi) 

Therefore  one  may  alternatively  compute  »’y(wi),  etc.  as  follows. 
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»'y(^l)  =  ^(*iK,(wi)  +  *  1*1  +  3*0  +  (-1)*0  =  1 

since  ty{xi)  =  1,  <y(x2)  =  3,  and  =  —1.  And  similarly,  for  i  =  2,  3  and  4, 

»'y(^2)  =  *y(*lK.(^2)  +  =  1*2  +  3*0  +  (-1)*0  =  2 

^yi^z)  =  ^(*l)*Xi('W3)  +  ^i^2)*4'^z)  +  tyi^zKS'^z)  =  1*1  +  3*0  +  {-1)*0  =  1 

*'y{^i)  =  +  ly(*2)«x,(w4)  +  ^(X3)«*,(W4)  =  1*0  +  3*1  +  {-l)*0  =  3 

Continuing  in  this  fashion,  we  obtain  ry(ws)  =  6,  ry(w8)  =  3,  ry(w7)  =-l,  ry(w8)  =-2,  and 
•■yKg)  =-l. 

In  order  to  determine  s  Q  t,  we  need  to  compute  for  i  =  1,...,9; 

»-y(wi)=  V  {<y(x)+«^(wi):Wi€S(x)} 
x€T(y) 

Equivalently,  ry{wj)  = 

{ly(*i)+«xi(’Wi):Wi€S(xi)}  V  {yxoj+s^Wj):  Wi€S(xo)}  V  {ly(x3)+s^^(Wi):  WjeSCxs)}. 

Now  Wj,  Wj,  W3  are  elements  of  S(xi),  but  none  is  an  element  of  either  ${x^  or  S{x3).  Thus 

•  ’•y('«^i)  *  <y(xi)+«xi(w)i)  =  1+1=2 

ry(w2)  =  «y(xi)+«,^(w)2)  =  1+2=3 

»'y('W3)  =  ly(Xi)+«x,(w)3)  =1+1=2 

Similarly,  since  w^,  W5,  Wg  €  5(x2)  and  Wy,  Wg,  W9  €  S(x3),  it  follows  that 

=3+1=4 

’■yl'^'s)  =  ly(X2)+«x,(w)5)  =3+2=5 
»’y('«'8)  =  ly{x2)+«xj('«')6)  =3+1=4 

ry(w7)  =  ty(x2)+«x,(w)7)  =  -1+1  =  0 
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''y('»^8)  =  ^(*3)+«x,('w)8)  =  -1+2  =  1 

*'y(’Wg)  =  ty(x3)+«,,(w)8)  =  -1+1  =0 

For  8®t  we  need  only  replace  the  +sign  in  the  definition  of  0  by  multiplication. 
Specifically, 

^(Wl)  =  ly(XiK/Wij  =  1*1  =  1 

ryM  =  <y(Xl)«Xj('W2)  =  1*2  =  2 
Similarly  for  ry(Wi)  with  i  =  3, ..,9. 

We  observe  that  in  this  particular  example,  we  have  the  coincidence  that  8@t  =t®8. 
The  reason  for  this  is  that  for  each  distinct  pair  x,x'€T(y),  S(x)riS(x')  =0.  If,  for  example, 
T(y)  were  a  von  Neumann  or  cruciform  configuration,  then  equality  need  not  hold. 

Subtraction,  division,  minimum,  scalar  multiplication  -  etc.,  can  be  defined  from  these 
basic  operations  in  a  straight  forward  manner.  Thus,  for  example,  t  •  s  is  defined  by  (t— = 

8y. 

Template  composition  and  decomposition  are  the  primary  reason  for  introducing  opera¬ 
tions  between  generalized  templates.  Composition  and  decomposition  of  templates  provide  a 
tool  for  algorithm  optimization.  For  instance,  if  s  and  t  are  as  in  the  example  above  and  r  = 
8  ©  t,  then  computation  of  a  @  r  =  a  @  (a  ©  t)  by  (a  ©  s)  ©  t  uses  6  local  multiplications 
instead  of  9. 

In  general,  if  r  is  an  nxn  template,  and  s  and  t  are  decompositions  of  r  into  IX  n  and 
nX  1  templates,  respectively,  then  the  computation  of  a  ©  r  by  (a  ©  s)  ©  t  uses  2n  multi¬ 
plications  instead  of  n^.  General  methods  for  template  decomposition  and  applications  of 
decompositions  to  algorithm  optimization  can  be  found  in  Reference  5  and  in  Section  III. 

10.  MULTIVALUE,  MULTIDATA  AND  MULTISENSOR  IMAGES 

n 

Let  a  be  an  F  valued  image  on  X.  If  X  is  homogeneous  and  F  =  IT  Fj  is  homogeneous 

i-l 

with  n  >  1,  then  a  is  called  a  multivalue  image.  If  F  is  heterogeneous,  then  a  is  called  a 
multidata  image. 
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EXAMPLE  2.14. 

A  LANDS-AT  image,  a,  of  n-spectral  bands  provides  an  example  of  a  multivalue 

a 

image.  Here  F  =  H  F;  =R"  and  a(x)=(aj(x), . ,an(x)),  with  each  a^(x)€R. 

n 

If  X  is  heterogeneous  and  F  =  H  Fj  is  homogeneous  with  n  >  1,  then  a  is  called  a 

i— 1 

multisensor /multivalue  image.  If  both  X  and  F  are  heterogeneous,  then  a  is  called  a 
multisensor/multidata  image.  If  X  is  heterogeneous  and  the  underlying  value  set  F  is 
known,  then  a  is  simply  called  a  multisensor  image. 

n 

If  F  =  n  Fj,  and  a  Ls  a  multivalue  or  multidata  image  on  X  with  values  in  F,  then  a 

can  always  be  viewed  as  a  stack  of  n  images  a=(ai,a2, . ,a„),  where  the  i-th  coordinate 

image  a^  of  the  stack  is  defined  as 

=  pi(a)  =  { (x,ai(x)) ;  ai(x) = Pi(a(x)) }. 

Multisensor  images  are  special  cases  of  F  valued  images  on  X  viewed  from  a  different 
perspective.  For  example,  let  X=XjX  X2,  where  Xj  and  Xo  denote  two  finite  rectangular 
arrays  in  shown  in  Example  2.4.  Now  suppose  that  aj  :  X,  — ►Fj  defines  an  Fj  valued 
image  on  Xj,  and  :  Xj-^Fj  defines  an  Fj  valued  image  on  X2;  i.e.  let 

.»!  =*  {(xi,ai(xi))  ;  XjGXi},  and  a,  *  {(X2,a2(x2))  :  X2€X2} 

Then  the  two  images,  aj  and  ao,  viewed  as  an  ordered  pair,  define  an  F  =  FjXFo  valued 
image  a  =  {&i,aLo)  on  X  =  XjXXn  if  we  define 

a:  X-^F  by  a(x)  =  a(x„X2)  =  (ai(xi),a2(x2)) 

That  is, 

a  =  {(x,a(x)):x€X}  =  {((xi,X2),(ai{xi),a.,(x2))  :  XiEXj,  X2€X2}. 


In  general,  if  i=l,2,...,n,  and  a^GF;  ,  that  is,  ai  =  {(Xj,aj(Xj)):XiGXi },  then  the  image  a 
defined  as 

a  =  {(x,a(x)):a(x)=((ai(xi),a2(x2)),  ■  •  ■  ,a„(x„))} 

n  n 

is  an  element  of  F^,  where  F  =  fl  F-,  and  X=  11  Xj.  Observe  that 

i— 1  i—l 


^  "  (nx,)  (nx,) 

F^=(nFi)^=(nFi)'-‘  =F'-‘  . 

i— 1  i— 1 
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However,  the  special  kind  of  multisensor  image  just  described  can  also  be  viewed  as  an  ele- 

ment  of  11  (Fj)"^  since  it  can  be  written  in  the  form  a=(ai,a2, . ,a„),  where  a^  €  Fj  .  We 

note  that  each  a^  may  again  be  a  multidata  or  multivalue  image  as  no  restriction  was 
imposed  on  the  value  set  Fj. 

11.  MULTIVALUE,  MULTK)ATA  AND  MULTISENSOR  IMAGE  OPERATIONS 

The  unary  and  binary  operations  on  multivalue/multidata  and  multisensor  images  are 
the  operations  induced  by  the  set  0  of  operations  of  the  underlying  value  set  F.  For 

instance,  if  F  =  11  F;,  a  =  {(x,a(x)) :  a(x)  =  (ai(x) . ,a„{x)) }  and  b  =  {{x,b(x))  : 

1—1 

b(x)=(bi(x), . ,b„(x))}  are  F  valued  images  on  X,  and  o=(oi, . ,o„)6  0  a  binary  opera¬ 

tion  on  F,  then 

a o  b  ={(x,c(x)) ;  c(x)  =(aj(x)  Oj  b,(x), . ,a„(x)  o„  b„(x)) } 

or,  equivalently, 

aob=(aiOibi, . ,a„  o„b„). 

If  o  is  unary,  then 

o(a)  =  {(x,c(x)) :  c(x)=(oi(ai{x)) . ,On(a„(x)))}. 

or,  equivalently, 

o(a)=<Oi(aJ, . o„(a„)). 

Example  2.15. 

If  F  =R"  and  o  is  a  binary  homogeneous  operation  with  o-,  =  +,  then 
a+b  =(a,  +b„ . a„  +b  J  =  {{x,c(x) ;  c(x)  =(a,(x)  -f  b,(x), . a„(x)  -t-b„(x)) } 

Similarly,  if  o  is  a  unary  homogeneous  operation  with  o-,  =  sin,  then 
sin(a)=(sin(ai) . ,sin(an))={(x,c(x);c(x)=(sin(ai{x)) . ,sin(a„(x)))}. 

and  if  Oj  =  x>)  then 

X>b(a)  =  {(x.c(x)) :  Ci(x)  =  1  if  ai(x)  >  bi(x),  otherwise  Ci(x)  =  0,  x  €  X} 
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If  F— R"  as  in  Example  2.15,  then  the  operations  defined  for  R  valued  images  extend  in 
a  natural  way  to  R"  valued  images.  For  example,  if  a  and  b  are  R"  valued  images  on  X, 
then  the  dot  product  is  defined  as 

a«b  =  S  a(x)b‘(x) 

x€X 

where  b^x)  denotes  the  transpose  of  b(x).  The  notion  of  scalar  multiplication  and  addition 
are  extended  to  vector  multiplication  and  vector  addition  in  a  similar  fashion.  In  particular,  if 
b={(x,b(x):  VxGX  be  viewed  as  a  vector 

b  =  k  =  (kj,  .  .  .  ,  kn)  €  R“  and  we  define 

ka  =  b*a  =  {(x,c(x)) :  Ci(x)  =  kiai(x),  x  €  X}. 
and 

k  +  a  =  b4a  =  {(x,c(x)) :  Cj(x)  =  ki4ai(x),  x  €  X  } 

Thus,  if  b  is  the  constant  vector  valued  image  k,  i.e,  b  =  k  =(k,...,k),  then  again, 

b*a  =  va  s  ka  and 
b4-ask+a  =  k+  a 

Obviously,  these  notions  'extend  to  F  =  C"  valued  images  and  vector  valued  images  in 
general.  In  comparison  to  the  next  set  of  operations,  the  operations  just  discussed  are  not 
new  operations,  but  simply  the  naturally  induced  operations  of  {F,0}. 

k  n 

Let  F  =  n  F;  and  F'=  IT  F':  be  two  value  sets.  If  k  ^  n,  then  any  function  g-F  ""^F^  is 

i-l  i-1 

called  a  data  splitting  function.  If  k  >  n,  then  g  is  called  a  data,  fusion  function.  In  either 
case,  g  is  of  the  form  g=(gi,g2,  •  •  •  >  gn)  where  gj  =  pjg. 

Example  2.16. 

The  function  g;  R  -♦  R^  is  defined  by  g(r)  =  (cos  r,  sin  r)  is  an  example  of  a  value 
splitting  function.  The  function  f:  R"  — *•  R  defined  by  f(ri,r2)  =  rj+To  provides  an 
example  of  a  data  fusion  function. 

One  of  the  most  common  data  fusion  operations  is  the  reduce  operation.  Given  a  homo- 
n 

geneous  value  set  F  =  n  F;,  a=(a,,  ....  a^)  an  F  valued  image,  and  o=(oi . ,On_,)  with 

i-l 

each  a,  a  binary  operation  on  Fj,  then  the  reduce  operation  /o(a)  is  defined  as 
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/o(a)=aiOi(ao02( . 


Example  2.17. 

Suppcee  a=(ai,  .  .  .  ,  a„)  an  F  valued  image  and  F=R“,  then 

k 

/"f{a)  =:5ja^  =  aj  +  Oo  +  ....  +  sl^ 
i— 1 

and 

/V(a)  =  Va,=ai  VaoV....  Va^ 

i— 1 

Observe  that  the  vector  valued  image  a  was  reduced  to  a  real  valued  image. 

Two  elementary  fusion/splitting  functions  of  prime  importance  are  the  projection  and 

n 

injection  functions.  Let  F  =  11  F:.  The  projection  onto  the  ith-coordinate  is  defined  as 

i-l 

before  by 

Pi  •  ^  ^i)  (*’l»  •  •  •  I  *’n)  • 

For  F  a  homogeneous  value  set  we  define  the  injection  function 

q:Fk-*-F,  by  q(r)=(r,r,...,r) 

The  ith-coordinate  injection  q-,  :  F  — ►  F  is  defined  as  qj  =  qp;.  Thus  qj  replaces  all  the 
coordinate  values  of  a  point  by  the  value  of  the  ith-coordinate,  namely 

qi(ri,r2,  .  .  .  ,rj,  .  .  .  ,r„)=(ri,ri,  .  .  .  ,  r;) 
n 

If  F  =  n  F:  and  F;  =  R  for  some  j  between  1  and  n,  then  F^  forms  a  lattice.  In  partic- 

i-l  ‘ 

ular,  we  define  the  jth-coordinate  maximum  and  minimum  of  two  F  valued  images  as 

i.  aVijb  =  a*qj[)^{a)l  -I-  b*qj[;)(^(a)l 
and 

ii.  aAljb  =  a*qj[x^(a)j  +  b*qj[x^(a)l, 
respectively,  where  [X>b(3-)1  denotes  the  Boolean  complement  of  )(>b(a).  Thus, 
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aVi jb  =  {(x,c(x)) :  c(x)  =  a(x)  if  aj(x)  >  bj(x),  otherwise  c(x)  =  b(x)  }. 

Ebcample  2.18.  Directional  Edge  Detection 

Directional  edge  detection  is  accomplished  by  convolving  the  image  with  the  follow¬ 
ing  six  3x3  edge  masks  with  each  mask  having  a  direction  associated  with  it. 
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Figure  13.  The  Six  Directional  Edge  Detection  Masks 

The  resulting  six  images  are  then  fused  to  form  a  single  2-valued  image  b.  Data 
fusion  is  accomplished  by  assigning  to  the  resultant  pixel  the  value  of  the  largest 
magnitude  of  the  corresponding  pixels  in  the  six  images  and  either  assigning  the 
direction  9  associated  with  the  mask  of  the  convolved  image  if  the  specific  pixel 
value  of  the  convolved  image  is  positive,  or  ^  -t-  180°  mod  360°  if  the  value  is 
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negative.  Thus,  each  pixel  of  b  has  both  a  magnitude  and  a  direction  associated 
with  it.  The  Image  Algebr?  translation  of  this  algorithm  is  as  follows.  Let  0-^ 
denote  the  direction  in  degrees  associated  with  m;.  Let  f:  R  — ►  be  defined  by 

f(r)  =(|rl,180Xc;o(r)),  then 

i.  aj  =  (0,^;)  +  f(a@mi),  i=l,...,6  and 

ii.  b=(Vli.)  a, 


The  next  three  figures  provide  a  pictorial  example  of  this  algorithm. 


Figure  14.  Input  Image  a 
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Figure  15.  Magnitude  Image  Pi(b) 
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h(2)  ~  ZjX  . X  Zj 

In  an  analogous  fashion,  let  (ai,a2, . ,an)^  n  (Fj^)  and  define  g:  11  (Fj^)  — ►  F^  by 

i— 1  i— 1 

g(ai»a2, ,an) =a={(x,a(x)) ;  a(x)  =(a,(xi),a2(x2), a„(xj) } 

&  n  Y 

Observe  that  g  is  a  one-to-one  and  onto  function  from  IT  (Fj  )  to  g(n  (Fj  ))CF^;  i.e., 

i«l  i— 1 

g  embeds  multisensor  images  into  F^.  Similarly,  h  is  a  one-to-one  and  onto  function  from 

n  (2^)  to  h(  n  (2'^)).  A  template  teFy^x  is  called  a  multilevel  template  if  its  range  is 
i— 1  i— 1 

h(n  (2^))  X  g(n  (Fj^))  c  2^x  F^ 

i— I  i— 1 

Y 

Given  a  multilevel  template  t  =  (T,  f),  we  can  define  T  :  Y  — ►  2"^  by  Tj  =  pj  h“^  T, 

n  V  V  Y 

where  pj  denotes  the  ith-coordinate  projection  IT  (2"^)  — ►  2^  and,  similarly,  ^  :  Y  — +■  F;  by 

n  Y  X. 

fc  =  p'ig"^  t,  where  p'j  denotes  the  ith-coordinate  projection  IT  (Fj  )  —*■  Fj  .  Thus,  a  mul- 

i-l 

tilevel  template  t  can  be  thought  of  as  a  stack  of  templates  t={ti,t2, . ,t„),  where  each  ith- 

coordinate  t;  =(Ti  ,  ^)  €  (Fi)Y_;Q. 

Operations  between  multivalued  or  multisensor  images  and  multilevel  templates  are 
natural  extensions  of  previously  defined  operations  between  images  and  templates.  Suppose 

n 

a€F^  and  tSFy-^w  multilevel  template  t=(tj,t2, . ,tj,  where  W=  FI  Wj  and  each 

i— 1 

pair  WjjXj  being  subsets  of  the  same  coordinate  space.  Let  o=(oi, . ,o„),  where  each  Oj  is 

an  operation  between  Fj  valued  images  on  and  Fj  valued  templates  from  Y  to  Wj.  Then 
aot  is  defined  as 

aot=(aiOiti, . ,a„  o„tJ. 

Ebcample  2.19. 

Suppose  F  =R^  and  o=(@,  Q  ,®  ),  then 

a  o  ^  —  (a^  ,  ao  0  to ,  a^  tgj 

If  o  is  homogeneous,  say  o  =(  ® ,  @ ,  © ),  then  -  following  our  earlier  convention  - 
we  set  o  =@  and  define 
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a@t=(ai@ti,a2@t2,ai@t3) 


It  follows  from  the  above  discussion  that  a  multilevel  template  can  have  different 
configurations  at  different  levels  and  can  operate  differently  on  the  different  levels  of  a 
multivalue  or  multisensor  image. 
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SECTION  IV 

RELATION  OF  IMAGE  ALGEBRA  TO  OTHER  MATHEMATICAL  STRUCTURES 


One  of  the  main  tools  of  algebra  is  the  isomorphism.  An  algebraic  structure  can  be 
thought  of  as  a  set  of  symbols  and  rules  for  combining  the  symbols.  The  structure  can  be 
used  to  model  systems  which  obey  those  rules.  If  two  algebraic  structures  are  isomorphic, 
that  is,  if  there  exists  a  1-1,  operation  preserving  mapping  of  one  onto  the  other,  then  they 
provide  two  different  viewpoints  of  the  same  situation.  The  idea  is  that  the  more  ways  one 
has  to  look  at  a  problem  the  better  chance  there  is  of  solving  it. 

Developing  an  algebraic  structure  for  digital  image  processing  enables  one  to  take 
advamtage  of  isomorphisms.  In  this  chapter  we  establish  isomorphisms  between  the  image 
algebra  and  the  mathematical  structures  known  as  linear  algebra,  polynomial  algebra,  and 
lattice  algebras.  We  also  show  how  these  relationships  can  be  applied  to  the  problem  of 
developing  systematic  techniques  for  the  optimization  and  derivation  of  parallel  algorithms. 

1.  RELATIONSHIP  BETWEEN  IMAGE  ALGEBRA  AND  LINEAR  ALGEBRA 

A  subalgebra  of  the  image  algebra  is  a  subcollection  of  the  set  of  operators  and  operands 
of  the  image  algebra.  In  this  section  a  precise  relationship  between  a  subalgebra  of  the  image 
algebra  and  linear  algebra  will  be  described  which  shows  how  all  the  operators  and  operands 
of  finite  dimensional  linear  algebra  (matrices,  vectors,  scalars  and  the  operations  between 
them)  are  modeled  by  the  image  algebra.  One  consequence  of  the  relationship  is  that  any 
linear  image- to- image  transform  can  be  written  as  a@k  for  some  template  k.  It  will  then  be 
shown  how  the  above  lelationship,  along  with  the  techniques  of  template  decomposition  and 
configuration  restriction,  can  be  used  to  develop  algorithms  for  implementing  linear 
transforms  of  images  on  arrays  of  processors  with  various  interconnection  schemes. 

Let  X  be  an  arbitrary  bounded  subset  of  Z*'  and  F  a  field.  Observe  that  if  s,t€FY_x 
such  that  Sy(x)  =  ty(x)  for  every  x€X  and  y€Y,  then  a©t  =  a©8  for  every  a€  Fy^x- 
Thus,  padding  a  configuration  with  zeros  has  no  effect  on  the  computation  of  a@t.  We  will 
call  two  templates  equivalent  if  the  above  is  true.  This  definition  of  equivalence  imposes  an 
equivalence  relation  on  Fy_x-  The  collection  of  equivalence  classes  forms  a  ring  (called  a 
quotient  ring  of  Fy_x)  which  we  will  denote  by  Ly_«.x-  will  call  Ly_*x  the  ring  of  tem¬ 
plates  with  minimal  configuration.  Practically  speaking,  Ly_x  =  {  t  G  Fy_x  •  ^(*)  #  0  if 
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and  only  if  x  €  T(y)  }  and  if  two  elements  of  Ly_»x  combined  using  +  or  ©  the  result  is 
again  in  Ly_»x>  minimal  configuration.  Denote  by  0  the  template  in  Ly_»x 

empty  configuration  at  every  point.  Then  0®t  =  t®0  =  0  for  every  template  t  €  Ly_»x- 
Furthermore,  for  s,  t  G  Ly_,x>  t  —  s  =  0  if  and  only  if  t  =  s  (This  property  does  not  hold 
for  Fy_*x)-  Whenever  we  consider  ®  we  will  assume  that  the  template  has  minimal 
configuration. 

Since  X  is  a  finite  set  it  can  be  linearly  ordered.  Thus  we  can  write  X  =  {xo,Xj,...,x„_i}. 
We  define  a  mapping  i/ :  F^— ►F®  by  u{a.)  =  (a(xo),a(xi),  ...  ,a(x„_i))‘.  Since  t{ra-+«b)  = 
ri^a)4si^b)  and  u  is  1-1  and  onto,  y  is  a  vector  space  isomorphism. 

Let  Lx=Lx_x  denote  the  ring  of  n  x  n  matrices  with  entries  from  F 

under  matrix  multiplication  and  addition.  For  any  t  G  we  define  a  matrix  Mj  =  (mjj) 
where  mjj  =  ix,(*i)-  Note  that  the  i^**  row  of  Mj  is  Define  a  mapping  4':Lx 

't'(t)=Mt. 

Theorem  4.1.1.  'I'  is  a  ring  isomorphism  of  (Lx  ,©,+)  onto  (A/„,*,-f-).  That  is,  if  t,8  G  Lx 
then 

i. )  'l'(8+t)  =  ^(8)-l-'I'(t)  or  M,^  =  , 

ii. )  ^(8®t)  =  4^(t)'I'(8)  or  =  MjM,. 


Proof:  1.)  Let  t,8  G  Lx-  Then 


'l'(8-rt)  = 


Since  i>((«-t-t),j|)  =  +  K^.)>  conclude  that  '1'(8-K)  = 

'l'(s)-l-'l'(t). 

2.)  Let  r=8@t.  Then,  by  definition  of  ®,  r^,(Xj)  =  S  t^(x)«^(xj).  By  definition  of 

'  '  ^  x€T(xi)  ■  ’  ^ 
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n— 1 

matrix  multiplication  M^M,  =  (cjj)  where  Cy  =  But  since  ^|(xk)  =  0 

if  Xjj  ^71(Xi)  we  must  have  Cy  =  »‘x,(xj)  which  implies  that  'li'(8@t)  = 

'if  is  clearly  onto.  Furthermore  is  1-1  since  if  'F{t)  =  '^(s)  then  'F(8— t)  =0.  But  0 
is  the  only  template  with  all  zero  grey  values  and  minimal  configuration.  Hence 
(s— t)  =  0  so  8  =t. 

Q.E.D. 

Observe  that  if  D  is  any  n  x  n  diagonal  matrix  then  t  =  'F“^(D)  has  the  property  that  if 

b  =  Ij  then  for  every  a  €  F^,  a*b  =  a©t.  Hence  diagonal  matrices  correspond  to  tem- 
x€X 

plates  which  are  essentially  images,  that  is,  the  image  b  can  be  used  to  implement  the 
transform  a  i— a@t. 

We  now  prove  a  final  theorem  concerning  the  relationship  between  image  and  linear 
algebra. 

Theorem  4.1.2.  For  every  a€  and  t  €  Lx,  a©t  =  i/~^('F(t)i^a)). 

Proof:  Let  b  —  a©t.  Then,  by  definition  of  ®,  b(y)  =  S  a(x)lL(x).  'F(t)i^a)  = 

x€T(y) 
n— 1 

(co,Ci,...,Cn_if  where  Cj  =  ^txi(X|t)a(Xfc).  Since  =0  if  x^  ^  71(Xi)  we  have  that 

Cj  =  b(xi)  which  shows  that  i^b)  =  ^(t)t{a).  The  conclusion  follows  from  the  fact 
that  1/  is  invertible. 

Q.E.D. 

The  above  theorems  have  many  consequences.  Two  important  ones  are  that  template 
composition  corresponds  to  matrix  multiplication  and  that  the  action  of  a  template  on  an 
image  via  the  ©  operation  corresponds  to  the  action  of  a  matrix  on  a  vector.  Thus  a 
subalgebra  of  the  image  algebra  serves  as  an  alternative  algebraic  formulation  of  finite 
dimensional  linear  operations.  By  incorporating  the  notion  of  configuration  in  the  definition 
of  templates  a  geometric  component  has  been  associated  with  a  linear  computation.  This 
concept  becomes  important  when  trying  to  map  algorithms  to  special  architectures  with  cer¬ 
tain  geometric  configurations  associated  with  them.  Since  these  theorems  are  true  for  any 
bounded  X  C  Z*"  all  linear  transforms  defined  on  any  of  the  examples  of  the  previous  section 
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can  be  expressed  using  ®.  Moreover,  all  decompositions  of  linear  transforms  can  be  expressed 
also. 

The  special  case  that  X  is  an  m  x  n  array  is  of  particular  interest.  We  take  X  =  {  {x,y) 
:0<x<m-l,0<y  <n-l  }.  We  can  order  X  lexicographically  in  row  major  form  by  map¬ 
ping  (x,y)  I— ►  xn-+-y.  The  mappings  u  and  'I'  then  have  the  form;  t{a)  =  (aQ,aj,...,ann,_i)'' 
where  a^  =  a(x,y)  if  i=xn-l-y  and 

m  = 

Henceforth,  whenever  X  is  an  m  x  n  array  we  shall  assume  that  these  particular  mappings 
have  been  used. 

2.  SIMULATION  OF  ARRAYS  OF  PROCESSORS 

The  image  algebra  can  be  used  to  simulate  an  array  of  processors  with  limited  commun¬ 
ication  by  imposing  restrictions  on  the  template  configurations  allowed  in  translating  algo¬ 
rithms.  In  fact,  suppose  an  array  of  processors  is  modeled  as  a  graph  G  =  (V(G),E(G))  where 
the  vertices  represent  the  processors  and  the  edges  represent  a  communications  link  between 
processors.  Suppose  G  is  embedded  in  n-dimensional  space.  Then  let  X  =  V{G)  and  define  a 
template  r  =  r(G)  €  Fx_x  by  ^(y)  =  {  x  :  [y,x)  €  E(G)  }  and  =  {(x,ry(x))  :  ry(x)=l  for 
every  x  €  ^(x)}.  We  say  that  an  expression  p(aj,ao,  •  •  ,a„)  is  compatible  with  the  tem¬ 
plate  restriction  induced  by  r  if  for  every  template  t  appearing  in  p  we  have  that 
T(y)  C  ^(y)  for  every  y  €  X  In  other  words,  if  p  is  compatible  with  r  then  the  algorithm 
defined  by  p  can,  in  principle,  be  implemented  on  the  processor  array  modeled  by  G.  In  this 
way  the  image  algebra  can  serve  as  an  algebraic  model  for  special  architectures. 

For  example,  assume  that  X  is  an  n  x  m  array  and  define  the  elementary  von  Neumann 
template  v  by  V(x,y)={(x+l,y),(x-l,y),(x,y-»-l),(x,y-l),(x,y)}  and  U(x,y)(u,v)  =  1  if  and  only  if 
(u,v)  G  F(x,y).  Then,  using  v  as  a  restriction  template,  the  image  algebra  simulates  a  mesh 
connected  array  of  processors.  Thus,  if  p(ai,a2)  =  a,©ti— 2(ao  (3  to)-+-b*(ai*ao®tj)  and 
Ti(x,y),72(x,y)  C  for  all  (x,y)  G  X  then,  in  principle,  the  algorithm  defined  by  pjaj.a^) 


53 


could  be  implemented  on  such  a  machine. 

Another  example  is  the  hierarchical  or  pyramidal  arrays.  Let  n  be  any  positive  integer 

n 

and  for  i=0,l,...,n  let  Xj  =  {  (x,y,z) :  1  <  x,y  <  S®"*  ,  z  =  i  }.  and  X  =  UXj.  Define  the  res- 

i-O 

triction  template  r  €  Fx_*x  =  {  (u,v,w)  : 


w  =  z  and  (u,v)  is  an  ^-neighbor  of  (x,y);  or 


w  =  z-l,x  = 


,  and  y  = 


or 


w  =  z+1,  u  =  1— 


,  and  V  = 


}nx. 


The  notation  fa]  denotes  the  smallest  integer  greater  than  or  equal  to  a.  An  image  algebra 
expression  p  over  X  will  be  compatible  with  this  pyramidal  structure  if  for  every  template  t 
appearing  in  p  we  have  T(x,y,z)  C  .2(x,y,z)  for  every  (x,y,z)  G  X  Thus,  if  n  =  2  then  the 
configuration  of  a  template  at  the  point  (1,1,1)  must  consist  of  a  subset  of  the  following 
points: 


the  8-neighbors  in  the  same  plane,  (1,2,1),  (2,2,1),  (2,1,1); 

the  children,  (3,3,0),  (3,2,0),  (3,1,0),  (2,3,0),  (2,2,0),  (2,1,0),  (1,3,0),  (1,2,0),  (1,1,0);  and 
the  parent,  (1,1,2). 


Any  image  algebra  expression  over  X  involving  only  templates  satisfying  the  template  res¬ 
triction  induced  by  r  defines  an  algorithm  which  could,  in  principle,  be  implemented  on  such 
a  machine. 


As  a  final  example  consider  the  binary  n-cube,  Reference  14.  It  can  be  represented  as 
the  set  of  n-tuples  with  coordinates  either  0  or  I,  that  is,  X  =  {  (xjjXo,...,  x„)  :  Xj  G  {0,1}  }. 
Define  the  restriction  template  r  by  declaring  r  to  be  the  elementary  template  with 
configuration  function  given  by  ^(y)  =  {  x  :  H(y,x)  <  1  }  where  H  denotes  the  Hamming 
distance. 


These  examples  show  how  the  image  algebra  can  be  used  to  simulate  special  architec¬ 
tures.  Since  the  interconnection  scheme  of  any  processor  array  can  be  represented  by  a 
graph,  virtually  any  processor  array  can  be  represented  by  the  image  algebra.  Moreover,  by 
our  completeness  theorems.  References  12  and  10,  not  only  can  the  arrays  themselves  be 
represented  but  the  transformations  defined  on  the  arrays  can  be  represented  also. 


54 


3.  TEMPLATE  DECOMPOSITIONS  AND  MATRIX  FACTORIZATIONS 

We  can  combine  the  theorems  on  the  relationships  between  image  and  linear  algebra 
with  the  use  of  the  image  algebra  as  an  algebraic  model  for  various  processor  arrays  to 
mathematically  formulate  the  problem  of  mapping  linear  transforms  to  specific  hardware.  If 
a  processor  array  is  modeled  as  a  graph  G,  r  is  the  restriction  template  associated  with  G 
and  t  is  a  template  compatible  with  r,  then  we  say  that  t  and  Mj  are  local  with  respect  to  r. 

A  template  r^triction  will  result  in  a  particular  matrix  structure.  For  example,  the  von 
Neumann  restriction  will  impose  the  condition  that  a  matrix  corresponding  to  an  allowable 
template  will  have  a  certain  banded  structure.  Many  techniques  for  decomposing  matrices 
have  been  developed.  If  t  is  a  template  and  if  a  local  decomposition  of  M^  can  be  found,  then 
the  mapping  can  be  used  to  construct  a  local  decomposition  of  t.  Moreover,  because  of 
the  algebraic  relations  of  associativity  and  distributivity,  the  local  decompositions  will  result 
in  local  algorithms.  That  is,  if  Mt  ^^Mo,  and  ='I'~*(Mo),  then  'P~*(Mt)  = 

t  =  t2@ti  and  a0t  =  a0(t2@tj)  =  (a@t.2)©tj.  Also  if  M^  =  M14M2,  then  =  t  = 

tj  + and  a@t  =  a©(ti  4-  to)  =  a@ti  +  a©t2.  Hence  we  can  decompose,  or  factor,  a  tem¬ 
plate  t  into  sums  and  products  of  local  templates  by  decomposing,  or  factoring,  the 
corresponding  matrix. 

Observe  that  factoring  into  products  implies  that  the  result  of  a  computation  can  be 
operated  on  at  the  next  step  whereas  decomposing  into  sums  implies  that  two  images  must 
be  computed  separately,  which  requires  twice  as  much  storage  space.  Since  all  the  opera¬ 
tions  are  preserved  under  the  mapping  such  a  decomposition  will  result  in  an  algorithm 
for  the  local  computation  of  a0t.  In  this  way  local  algorithms  can  be  developed  using  the 
interaction  between  the  image  algebra  and  linear  algebra.  The  algorithms  can  be  thought  of 
as  being  mapped  onto  a  specific  hardware  using  4'“'.  Therefore,  we  can  use  the  image  alge¬ 
bra  to  mathematically  formulate  the  problem  of  mapping  linear  algorithms  to  specific 
hardware  as  follows: 

Given  a  template  t,  find  a  decomposition  of  t  in  terms  of  sums  and  products  of  local 
templates,  or,  equivalently,  given  a  matrix  M,  find  a  decomposition  of  M  in  terms  of  sums 
and  products  of  local  matrices. 

We  have  used  these  ideas  to  develop  methods  for  deriving  algorithms  for  the  computa¬ 
tion  of  two-dimensional  discrete  Fourier  transforms  of  arbitrary  size  on  mesh  connected 
arrays.  We  show  how  to  decompose  the  two-dimensional  Fourier  template  into  products  of 
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local  templates.  These  methods  will  be  discussed  in  the  subsequent  sections. 

4.  IMAGE  ALGEBRA  AND  POLYNOMIAL  ALGEBRA 

In  this  section  precise  relationships  between  polynomial  algebra  and  the  image  algebra 
atre  described.  It  is  shown  how  these  relationships  can  be  used  to  develop  systematic  methods 
for  designing  paradlel  algorithms  for  circulant  image  transforms.  The  methodology  is  similar 
to  that  of  the  previous  section  in  that  we  seek  to  decompose  templates  by  decomposing  other 
algebraic  objects. 

Throughout  this  section  it  will  be  assumed  that  X  is  an  m  x  n  array.  Everything  that 
will  be  done  generalizes  easily  to  the  case  that  X  is  a  rectangular  parallelepiped  but  the  nota¬ 
tion  becomes  cumbersome.  We  will  show  that  the  ring  of  circulant  templates  on  X  is  iso¬ 
morphic  to  a  quotient  of  the  ring  of  polynomiab  in  two  variables.  We  exhibit  a  class  of  map¬ 
pings  which  are  related  to  this  isomorphism  and  which  are  almost  isomorphisms  themselves. 

We  begin  with  some  preliminary  definitions  and  notation.  If  x  E  Z*^  then  x 
mod(ni,n2,.--,  =  (^i  modiii.Xo  modn2,...,Xij  modn^). 

Definition  4.4.1.  If  ^  :  X  X  is  of  the  form  ^x)  ==  (x-f^)  mod(m,n)  Then  <j>  is  called  a 
circulant  translation. 

Definition  4.4.2.  If  t  €  Fx_»x  satisfies  the  following  two  properties  for  every  circulant 
translation 

every  y  E  X  and 

*0(y)(<^x))  =  ^(x)  for  every  x,y  E  X 
then  t  is  called  a  circulant  template. 

A  circulant  template  wraps  around  at  the  boundary.  We  denote  the  ring,  under  the 
operations  @  and  -f,  of  all  circulant  templates  in  L^  by  C^.  If  t  is  a  circulant  template  then 
we  call  a  mapping  f  ;  of  the  form  f(a)  =  a©t  a  circulant  transform.  Observe  that 

a@t  denotes  the  circular  convolution  of  the  images  a  and  b  where  b(i,j)  =  t(o,o)(('‘''j) 
mod(m,n)).  Thus,  circulant  templates  are  used  to  implement  circular  convolutions  which  are 
standard  techniques  in  digital  signal  and  image  processing. 

Let  K[x,y]/(x'"-l,y''-l)  denote  the  ring  of  polynomials  in  two  variables  where  addition  is 
performed  in  the  usual  way  and  multiplication  is  performed  as  usual  except  that  the  substitu¬ 
tions  x"*  =  1  and  y“  =  1  can  be  made  at  any  time.  We  say  that  a  polynomial  is  in  reduced 
form  in  this  ring  if  all  possible  substitutions  of  that  type  have  been  made.  Two  polynomials 
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in  this  ring  are  declared  to  be  equivalent  if  they  have  the  same  reduced  form.  For  example 
let  m  s:  3  and  n  =  2.  Then  =  1  and  =  1.  Let  f(x,y)  =  1+x+yx"  and  g(x,y)  =  l+yx+x^. 
Then  f(x,y)  +g(x,y)  ==  2+x4yx+yx^+x^  =2+(14-y)x+(l+y)x^.  Also  f(x,y)g(x,y)  = 
l+yx+x+x^+x+yx“+x^4-yx^+y^+yx^.  The  reduced  form  of  f(x,y)g(x,y)  is 
l+yx+x^+x+yx^+l+yx'+l+yx  =  34{2y+l)x4{2y+l)x"  . 

For  every  z  €  X  we  define  a  mapping  F^  :  — ►  K.(x,y|/(x'"— l,y"— 1)  by 

m— In— 1 

r,(t)  S  Pt(x,y,z)  =  E  E  i,(i,j)xy. 

1^  j-0 

If  t  is  a  circulant  template  then  F,  is  called  a  polynomial  representative  of  t.  It  is  easy  to  see 
that  a  circulant  template  t  is  local  with  respect  to  a  mesh  connected  array  if  Pf(x,y,0)  =  x-a 
or  y-a  for  some  a.  We  also  point  out  that  if  the  polynomial  corresponding  to  a  circulant  tem¬ 
plate  t  is  of  the  form  x‘y^  then  the  circulant  transform  a  i— ►  a@t  simply  circularly  shifts  all 
the  grey  levels  i  units  vertically  and  j  units  horizontally. 

Theorem  4.4.3.  If  s  and  t  are  circulant  templates,  then  Fo(8@t)  =  FQ(8)FQ(t). 

Proof:  For  any  z  €  X  denote  by  the  circulant  translation  defined  by  ^^(x)  =«  (x-z) 
mod(m,n).  Note  that  <Aj(z)  =  0.  Recall  that,  if  r  =»  s@t,  then 

*■0  =  { (y.»-o(y)) :  »’o(y)  = 

z€T(0) 

Since  s  is  a  circulant  «j(y)  =  «o{^z{y))  ~  *o((y*®)  mod(m,n)).  But  by  definition  of 
polynomial  multiplication  (with  x"'=y"=l)  the  numbers  E  to(*)V{y~®)  inod(m,n)) 

167(0) 

are  precisely  the  coefficients  of  the  polynomial  product  Pg(x,y,0)pt(x,y,0). 

Q.E.D. 

Theorem  4.4.4.  For  every  z  =(j,k)  €  X,  F^  is  1-1  and  onto.  Moreover,  if  s.t  6  6^,  then 

F,(s)-l-F,(t)  =  F,(s+t),  (7) 
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r.(t)-*'y'‘ro(t),and 

(8) 

r.(ii®t)xiy'‘=r.(a)r.(t). 

(9) 

Proof:  Fj  is  clearly  onto.  F^  is  1-1  since  contains  only  circulant  templates  with  minimal 
configuration.  Equation  (7)  follows  from  the  fact  that  addition-  of  templates  is 
defined  pointwise. 

To  prove  (8)  let  <t>  be  the  circulant  translation  defined  by  <^x)  =  (x-l-*)  mod(m,n). 
Note  that  4{Q)  —  z.  Let  t  be  an  arbitrary  circulant  template.  Then 

m— In— 1 

xjy‘'Pt(x,y,0)  =  xjx*'  E  S  io(j.b)x'y‘’  = 

i-O  h-0 


m-ln— 1 

=  E  E  ij((i-H)  mod  m,(h-fk)  mod  n)x'''’V*'''’'‘  = 

M  h-O 


m— In— l 

=  E  S  ti(i,h)x‘y'’  =  Pt(x,y,z). 

hM) 


since  x'"  =  y"  =  1  and  t  is  circulant.  This  enables  us  to  deduce  that  F^(t)  ==  xjy*'Fo(t) 
or  x'"'jy“-‘‘r,(t)  =  Fo(t). 


To  prove  (9)  observe  that  since  FQ(8@t)  =  FQ(3)FQ(t)  we  have  that 
r,{80t)  =  xjy‘'Fo{8)Fo(t)  =  xV'  [x'"-V"-''F,(8)x-^y'’-''F,(t)]  =  x'"-jyn-kr^(3)rjt). 


The  desired  result  now  follows  by  multiplying  the  last  equation  by  x^y''  . 

Q.E.D. 

Corollary  4.4.5.  Fg  is  an  isomorphism. 

Corollary  4.4.8.  Let  Zj  =  (i,j),  z^  =  (s.t)  G  X  and  t  €  C^.  If  Pt(x,y.Zi)  =  p,(x,y)p2(x,y), 
then  Pt(x,y,Z2)  =  q,(x,y)q2(x,y)  where  qi(x.y)  =  x=“'y’'~Jpi(x,y)  and  q.2{x,y)  =p.2(xy). 
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Corollary  4.4.7,  If  *  =  (i,j)  €  X,  t  €  C^,  and  Pt(x,y,z)  =  Pi{x,y)p2(x,y),  then  t  = 
ti®t2©8  where  s  represents  a  circular  shift  of  n-j  units  horizontally  and  m-i  units  vertically. 
Moreover,  since  s  represents  a  circular  shift  it  can  be  replaced  by  a  template  s  which 
represents  a  circular  shift  j  units  horizontally  and  i  vertically. 

The  Fj ’s  constitute  the  class  of  mappings  between  image  and  polynomial  algebra  which 
are  almost  isomorphisms.  Since  they  differ  from  isomorphisms  only  by  shifts  they  can  be  used 
in  the  same  way.  Clearly  if  any  of  the  polynomial  representatives  of  a  template  can  be  fac¬ 
tored  the  template  can  be  factored  correspondingly.  We  now  examine  some  specific  conse¬ 
quences  of  these  theorems.  We  will  use  them  to  show  how  any  separable  circulant  template 
can  be  computed  locally  with  respect  to  the  von  Neumann  configuration  and  give  upper 
bounds  on  the  number  of  parallel  steps  required  to  do  so.  Since  the  von  Neumann  restriction 
on  an  m  X  n  array  simulates  mesh  connected  arrays  these  methods  can  be  used  on  such 
machines. 

a.  Parallel  Algorithms  via  Polynomial  Factorization 


m— In-l 

Definition  4.4.8.  If  f(x,y)  =  X!  I]  ^ijX'y^  then 

iaiO  j-0  ' 

deg(f(x,y))  =  max  {  i-f-j  :  0  }. 

degx(f{x,y))  =  max  {  i  :  a,j  0  }. 
degy(f(x,y))  =  max  {  j  ;  5^  0  }. 

Definition  4.4.9.  If  t  is  a  circulant  template  and  deg(pt(x,y,z))  <  deg(pt(x,y,w))  for  every 
w  €  X  then  we  say  that  z  is  a  minimal  point  for  t. 

Definition  4.4.10.  Let  t  be  a  circulant  template.  We  say  that  t  is  separable  if  there  exists 
polynomials  f(x)  and  g(y)  such  that  Pt(x,y,0)  =  f(x)g(y).  By  the  preceding  corollaries  if  t  is 
separable  then  for  every  z  €  X  there  exists  polynomials  f(x)  and  g(y)  such  that  pj.{x,y,z)  = 

f(”)g(y)- 

We  will  make  use  of  the  fundamental  theorem  of  algebra.  We  state  it  here  for 
reference. 
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Fundamental  Theorem  of  Algebra:  If  r(x)  is  a  polynomial  in  one  variable  with  com¬ 
plex  coefficients  then  f(x)  can  be  factored  into  n  linear  factors.  That  is,  f(x)  = 
a(x— TjXx— Tj)  •  •  •  where  a,rj,r2,  '  •  ,r„  are  (possibly)  complex  numbers. 

The  following  theorem  is  a  very  satisfying  result  which  yields  a  systematic  method 
for  computing  any  separable  circulant  transform  locally  with  respect  to  a  4-connected  proces¬ 
sor  array  and  gives  an  upper  bound  on  the  number  of  parallel  steps  required  to  do  so. 

Theorem  4.4.11.  Let  t  be  a  separable,  circulant  template  and  let  z  =  (s,t)  be  a  minimal 
point  for  t.  Let  k  =  degx(pt(x,y,z))  and  j  =  degy(pt(x,y,z)).  Let  v  =  min(s,m-s)  and  h  = 
min(t,n-t).  Then  the  circulant  transform  a  — ►  a@t  can  be  computed  in  at  most  v-f-h-l-k-f-j-f-1 
local,  parallel  steps.  v-Hi  of  these  steps  consist  of  vertical  or  horizontal  circular  shifts  of  the 
entire  array  by  one  location,  kH-j  of  the  steps  consist  of  at  most  one  addition  and  one  multi¬ 
plication  (pos-^ibly  complex)  per  pixel,  and  one  of  the  steps  consists  of  at  most  one  multiplica¬ 
tion  per  pixel. 

Proof:  Since  t  is  separable  there  exists  f(x)  and  g{y)  such  that  Pt(x,y,z)  =  f(x)g(y).  By 

assumption  deg(f(x))  =  k  and  deg(g(y))  =  j.  By  the  fundamental  theorem  of  algebra 
f(x)  =  a(x-qi)(x-q2)...(x-qk)  and  g(y)  =  b(y-ri)(y-r2)...(y-rj)  where 
a,b,qi,...,q|j,ri,...,rj  are  (possibly)  complex  numbers.  Choose  circulant  templates  q,  r, 
qi  ,-,qk-  ri,...,rj  such  that  ro(q)  =f(x),  ro(r)  =g(y),  ro(qi)  =x-qj,  and  ro(ri)  =y-ri. 
Then  rj(t)  =  rg(q)rQ(r).  By  Corollary  4.4.7,  t  =  q©r@8  where  s  represents  a  shift 

of  V  locations  in  the  vertical  direction  and  h  in  the  horizontal  direction.  Now  by 

k  k  j 

Theorem  4.4.3,  ro(q)  =  a(nr'o(qi))  implies  that  q  =  a(@q).  Similarly  r  =  b(@ri). 

i-i  ‘-i  '-1 

k  j 

Hence  t  =  ab  (@q)@(@rj)@s  .  Since  templates  corresponding  to  linear  polynomi- 
i— l  i— l 

als  in  one  variable  are  local  and  shifts  are  obviously  local  this  last  equation  expresses 
t  as  a  product  of  local  templates.  This  concludes  the  proof. 

Q.E.D. 

The  templates  are  shown  in  Figure  17.  In  this  case,  one  thinks  of  the  templates  as 

masks.  Thus,  convolving  an  input  image  with  the  masks  qi,q2,  .  .  .  ,q|f  will  result  in  the 

k 

same  output  image  as  convolving  the  input  image  with  the  k  x  1  mask  q  = .®  ^ 


qj  = 


~<li 

1 


Figure  17.  Templates  Required  to  Compute  a  Separable  Circulant  Transform  Locally 

As  an  example  suppose  X  is  a  512  x  512  array  and  t  is  a  30  x  30  separable  tem¬ 
plate.  Then  z  =(14,14)  and  v  =  h  =  14.  Hence,  if  a  6  F^,  then  a  — ►  a@t  can  be  computed 
locally  in  parallel  with  60  parallel  steps  consisting  of  at  most  one  multiplication  and  addition 
per  pixel,  1  step  consisting  of  one  multiplication  per  point,  and  a  total  of  28  unit  horizontal 
or  vertical  circular  shifts.  Since  the  computation  in  its  original  form  required  900  multiplica¬ 
tions,  it  is  clear  that  template  decompositions  can  be  used  to  derive  algorithms  which  are 
more  efficient  with  respect  to  the  number  of  arithmetic  operations  as  well  as  parallel. 

Nontrivial  examples  of  such  templates  are  the  discretizations  of  the  Marr-Hildreth 
edge  operators,  Reference  15.  Other  examples  include  a  class  of  templates  obtained  by  com¬ 
puting  the  inverses  of  blurring  transformations,  a  specific  example  of  which  will  be  presented 
in  the  next  section. 

If  a  circulant  template  is  not  separable,  then  things  do  not  work  out  as  nicely.  This 
is  because  the  fundamental  theorem  of  algebra  is  not  true  for  multivariable  polynomials. 
However,  the  following  unique  factorization  theorem  of  polynomial  algebra  applies: 

Unique  Factorization.  If  f(x,y)  is  a  polynomial  of  two  variables  with  coefficients  from  F, 
k 

then  f(x,y)  =  nPi(x,y)  where  the  Pi(x,y)  are  irreducible  over  F. 
i-l 

Hence,  if  X  is  an  m  x  n  array  and  t  G  C^,  then  for  any  z  €  X,  Pj(x,y,z)  can  be 

k 

thought  of  as  an  element  of  F[x,yj  and  can  therefore  be  factored  as  Pt(x,y,z)  =nPi(x,y). 

i-l 

Since  Pj  has  no  power  higher  than  m  or  n  the  factorization  remains  valid  after  making  the 

k 

substitutions  x‘"=l  and  y''=l.  Thus  t  =tj©to©  •  •  t|j©3  =(©tj)©s  where  the  tj  are  the 

i-l 

templates  corresponding  to  the  p;  and  s  is  a  circular  shift  template. 

Thus  a  nonseparable  circulant  template  can  be  decomposed  into  a  product  of  tem¬ 
plates  which  are  irreducible  with  respect  to  polynomial  factorization.  These  irreducible  fac¬ 
tors  can  then  possibly  be  handled  by  techniques  such  as  matrix  decompositions  and/or 
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FFT’s.  Algorithms  for  factoring  multivariable  polynomials  as  well  as  a  large  number  of  refer¬ 
ences  on  the  subject  can  be  found  in  References  16  and  17 

Since  circulant  templates  compute  circular  convolutions,  it  is  not  surprising  that 
polynomial  algebra  is  helpful  in  designing  pau'ailei  algorithms  to  compute  them.  Indeed,  much 
of  the  theory  of  fast  transforms,  particularly  the  work  of  Winograd,  References  18  and  19, 
and  in  Nussbaumer,  References  20  and  19,  is  based  on  the  use  of  the  Chinese  Remainder 
Theorem  for  polynomials  to  derive  systematic  methods  for  developing  fast  convolution  algo¬ 
rithms.  Our  approach  is  fundamentally  different.  Whereas  their  criteria  is  the  reduction  of 
the  number  of  arithmetic  operations,  particularly  multiplications,  our  criteria  is  the  decom¬ 
position  of  the  templates  required  to  execute  a  given  transformation  to  a  form  compatible 
with  a  certain  geometric  configuration.  We  use  the  tool  of  polynomial  factorization  rather 
than  the  Chinese  Remainder  Theorem. 

We  have  seen  that  these  criteria  can  be  compatible.  Specifically,  template  decompo¬ 
sition  usually  does  result  in  a  reduced  number  of  arithmetic  operations.  Thus  the  methodol¬ 
ogy  of  parallel  algorithms  via  template  decomposition  and  configuration  restrictions  often 
yields  arithmetically  efficient  parallel  algorithms. 

5.  PARALLEL  ALGORITHMS  IN  THE  MAGE'  ALGEBRA 

In  this  section  specific  applications  of  the  general  techniques  described  in  the  previous 
section  are  presented.  The  first  is  a  class  of  algorithms  for  computing  two-dimensional  DFT’s 
of  arbitrary  size  on  mesh  connected  arrays.  The  development  of  these  algorithms  is  based  on 
the  use  of  the  technique  of  template  decomposition  via  matrix  factorization.  The  second 
application  is  an  algorithm  for  computing  a  global  separable  circulant  template  locally  (with 
respect  to  mesh  connected  arrays)  which  uses  Theorem  4.4.11.  The  template  arises  as  the 
inverse  of  a  separable  blurring  template.  Thus  this  technique  can  be  used  to  develop  parallel 
algorithms  for  image  restoration  via  inverse  filtering  when  the  degradation  model  is  separ¬ 
able. 

a.  Two-Dimensional  FFT’s  of  Arbitrary  Size  for  Mesh  Connected  Arrays 

In  this  section  we  show  how  the  relationship  between  the  image  algebra  and  linear 
algebra  can  be  used  to  develop  local,  parallel  algorithms  for  the  two-dimensional  FFT  for 
images  of  arbitrary  size.  We  will  define  the  template  used  to  compute  the  discrete  Fourier 
transform  (DFT).  Since  the  DFT  is  separable  the  matrix  corresponding  to  this  template  can 
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be  factored  into  a  Kronecker  product  of  one  dimensional  Fourier  matrices.  These  matrices 
can  be  factored  into  a  product  of  tridia^onal  matrices  and  permutation  matrices.  The  permu¬ 
tations  can  be  done  in  parallel  by  the  odd-even  transposition  sort  while  the  tridiagonal 
matrices  are  in  fact  local  with  respect  to  mesh  connected  arrays.  Hence,  translating  the 
decompositions  back  into  the  image  algebra  will  result  in  parallel  algorithms.  We  illustrate 
this  procedure  for  the  case  that  X  is  a  100  x  100  array.  Pease  used  matrix  algebra  to  show 
how  radix  two  FFT’s  could  be  implemented  in  parallel,  Reference  21.  Radix  two  FFT’s  for 
mesh  connected  arrays  have  also  been  developed.  References  22  and  23. 

b.  Notation  and  Basic  Properties 

Unless  otherwise  stated  we  consider  matrices  and  images  as  ordered  from  0  to  n-1. 
We  denote  by  the  group  of  permutations  of  n  objects,  and  by  Z„  the  ring  of  integers 
modulo  n.  We  think  of  as  acting  on  Z„.  Let  m  €  Z  with  m  >  1  and  let  =  exp(-27ri/m) 
where  i  =(— 1)*/^. 


Definition  4.5.1.  Let  <7  €  and  define  the  n  x  n  matrix  P^  by 


P<r  =*  (Pij)  where  pjj  = 


1  if  j  =  <7(0 

0  otherwise 


P,,  is  called  a  permutation  matrix. 

Assume  for  the  rest  of  the  section  that  n  =  nin2. 

Definition  4.5.2.  Define  (T^,  :  Z^— ^Z^  by  the  following  rule: 

If  i  €  Zn  and  i  =  a-l-bni  with  0  <  a  <  Hj  and  0  <  b  <  n^,  then  (7„^„Ji)  =  ano-f-b. 
*^n,n2  ^  called  a  generalized  shuffle  permutation. 

Definition  4.5.3.  Define  P(n,,n2)  €  AL  by  P(ni,no)  =  .  P(ni,no)  is  called  a  generalized 

“i“s 

shuffle  permutation  matrix. 

Observe  that  P{nj,n2)  ==  P^(n2,ni)  =  P“*(n2,ni). 
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Definition  4.5.4.  Let  A  be  an  n  x  m  matrix  and  B  an  s  x  t  matrix  both  with  entries  in  F. 
Then  the  Kronecker,  or  Tensor  Product  of  A  and  B,  denoted  A0B,  is  the  nsX  mt  niatrix 
given  by 


A(8lB  = 


^-1,0^ 


It  is  well  known  that  (A®B)(C®D)=AC®BD  provided  that  the  matrices  are  ail  of 

'  k  1  k 


the  appropriate  dimensions.  Hence 


nAi0ij=n  j^Aj^Ijl  The  Kronecker  Product  is  not 
i-o  ^  ’ 


commutative.  Rose  has  shown  the  following  (Reference  24  ); 


Eacli  If  A  is  an  n  X  m  matrix  and  B  is  a  p  x  q  matrix  then  P(n,p)(A(8B)P(q,m)=B®  A.  In 
particular,  if  A  is  an  n  x  n  matrix  and  B  is  an  m  x  m  matrix,  then 
P(n,mXA®B)P(m,n)=B®A. 

Definition  4.S.S.  Define  the  n  x  n  matrix  D(ni,n2)  by 

D(ni,n2)  =  blockdiag[D,5(ni,n2)l,  k=0,l,...,n2— 1. 

where 

Dk(“i.n2)  =diag(w‘'‘),  i=0,l,...,ni-l,  k=0,l,...,n2-l 


For  any  p,k  >  1  we  denote  =  D(p'‘  *,p)  and  P^k  =  P(p'‘~^p). 

Definition  4.5.6.  Define  the  n  x  n  matrix  by  F„  =  F^  is  called  the  one¬ 

dimensional  Fourier  matrix  of  order  n.  If  x  €  F"  then  F„x  is  the  one-dimensional  discrete 
Fourier  transform  (1-dim  DFT)  of  x. 

Definition  4.5.7.  Let  X  be  an  m  x  n  array  .  Define  g  G  Fx_*x  by  ,$(u,v)  =X  and 
ff(u,y)  =  (nm)“'/^exp[— 27ri(uj/n-t-vk/m)],  where  j,k  €  F^  are  defined  by  j  =  {  (x,y),j(x,y))  : 
j(x,y)  =  X  }  and  k  =  {  (x,y),k(x,y))  :  k(x,y)  =y  }.  g  is  the  two  dimensional  DFT  template. 
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The  transform  a  i— ►  a@g  is  the  DFT  of  a. 

Definition  4.5.8.  Let  X  be  an  m  x  n  array  and  let  g  be  as  in  the  previous  definition.  Then 
let  =  'F(g).  F,^  is  called  the  two-dimensional  Fourier  matrix  of  order  m  x  n. 

Notice  that  in  general.  Notice  also  that 

Fmxn  =  F^0Fn  =  (F^® In)(Im® ^n)-  observation  is  an  equation  which  expresses  the 

fact  that  the  two  dimensional  DFT  can  be  computed  by  first  computing  one  dimensional 
DFT’s  along  each  row  and  then  along  each  column  of  the  result.  A  tridiagonal  matrix  com¬ 
putes  a  linear  transformation  of  a  linear  array  locally.  Permutations  of  linear  arrays  can  be 
computed  locally  using  the  odd-even  transposition  sort,  Reference  25.  Hence  if  tridiagonal 
decompositions  of  the  one  dimensional  Fourier  matrix  can  be  found  for  every  n  then  local 
DFT’s  could  be  implemented  on  mesh  connected  arrays.  These  decompositions  have  been 
derived.  We  present  them  in  the  next  section. 

c.  Tridiagonal  Decompositions  of  the  Fourier  Matrices 

In  this  section  we  present  the  sequence  of  theorems  resulting  in  the  method  for  fac¬ 
toring  the  Fourier  matrices  into  products  of  tridiagonal  matrices  and  permutation  matrices. 
The  proofs  can  be  found  in  Reference  6.  It  is  our  intention  here  to  show  how  they  can  be 
used  along  with  the  image  algebra  to  derive  algorithms  for  the  local  implementation  of 
FFT’s  on  mesh  connected  arrays.  The  statements  of  the  theorems  show  how  to  proceed  and 
therefore  we  present  them. 

Assume  that  n  =  mk.  Rose,  Reference  24,  has  shown  that 
Fn  =(F,0Ik)D(k,m)(I,0F,)P(k,m) 

and,  that  if  m  and  k  are  relatively  prime,  then  there  exists  permutation  matrices  Qj  and  Q2 
such  that 

F„  =Qi{F„,®FJQ2 

Using  these  results  one  can  prove  the  following  theorems: 

3  k 

Theorem  4.5.9.  Suppose  n  =  JJp; '  with  s  >  2.  Let 

i-l 
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i-l 

1 


if  1  <  j  <  s 

if  j  =  -1,0  ^  "/“j  ^  i  ^ 


and  qj  =  Pj  ^  for  1  <  j  <  s,  qo  =  1.  Then  there  exists  permutation  matrices  Qq,  Q8_2, 

H  such  that 


F 


a 


U-O  J 


Theorem  4.5.10.  F 

Gp 


pk  =  Gp(Ipk-i<8Fp)Hp ,  where 


and 


~  n  Kk-m®Ppni]- 

m-2  J 


The  previous  theorems  swe  matrix  identities  associated  with  FFT’s.  When  the  blocklength  is 
prime  the  FFT  no  longer  works.  This  is  one  reason  why  the  radix  two  case  is  much  easier  to 
develop  a  parallel  algorithm  for.  Indeed,  the  last  theorem,  with  p  =  2,  is  a  tridiagonal 
decomposition  of  power  of  two  Fourier  matrices  and  therefore  yields  a  local  algorithm.  For 
the  general  case  we  resort  to  other  matrix  identities.  These  identities  are  based  on  the  Rader 
prime  algorithm,  Reference  26,  and  the  circular  convolution  theorem  in  matrix  form. 

Some  notation  is  required.  Let  A  be  an  n  x  n  matrix.  We  denote  by  the  (n-Hn) 
X  (n-kn)  matrix 


Im  o' 
0  A 


We  denote  by  Up  the  p  x  p  matrix  with  all  one’s  in  the  first  column  and  one’s  down  the  diag¬ 
onal.  For  example  if  p  =  5  then 
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Us 


1  0  0  0  0 
110  0  0 
10  10  0 
10  0  10 
1  0  0  0  1 


Let  p  be  a  prime.  Then  it  is  well  known  that  for  every  positive  integer  a  with  1  <  a  <  p, 
p-1  is  the  smallest  positive  integer  such  that  (mod  p)  =  1.  Choose  any  such  a  and  let  cj 
=  1.  Define  the  polynomial  fp(x)  =  Co+CiX+  •  •  ■  -fCp.jx'*'^.  Then  define  the  p  x  p 

matrix  Ap  by  Ap  =diag(l,fp(a;p_i‘-i)),  i=l,2 . p-1. 

Theorem  4.5.11.  (Prime  Decomposition  Theorem)  Suppose  p  is  a  prime.  Then  there  exists 
permutation  matrices  and  Ro  such  that 

where  *  denotes  the  conjugate  transpose. 

Observe  that  if  A  and  B  are  any  n  x  n  matrices  then  {IjSiA*)  =  (I^0A)*  and 
=b(AB)^"*1  Therefore  if  “XITj  where  the  Tj’s  are  tridiagonal  matrices  then 

~  is  a  product  of  tridiagonals.  Thus 

the  tridiagonalization  problem  is  reduced  to  factoring  matrices  of  the  form  Up  and  Fp_/^^. 
But  Fp_i(^^  =  blockdiag[l,Fp_il  SO  it  IS  clcsir  thdrt  one  C3>ii  proceed  inductively  until  p—2  or 
p=3  in  order  to  factor  Fp_/‘l  It  is  easy  to  see  how  to  factor  Up  by  considering  an  example. 

Suppose  that  p  =  -5.  Then  if  x  =  (xo,Xi,X2,X3,x.,)*’,  U5X  = 
(xQ,Xo-t-Xi,XQ+X2,XQ-Px3,XQ-i-xJ.  The  algorithm  depicted  in  Figure  18  will  compute  UpX  locally 
in  p-1  steps.  In  matrix  form 
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0  Oi 


1  0 
1  1  0 


*0 
*1 
*2 
*3 
*4 

Hence  the  problem  is  successively  reduced  from  F„  to  FpH  to  Fp  to  Fp.^.  This  reduc¬ 
tion  process  terminates  when  p=2  or  p==3  in  which  case  we  have  all  tridiagonal  or  permuta¬ 
tion  factors.  We  now  show  how  these  theorems  can  be  used. 

d.  Permutations 

In  this  section  we  describe  how  separable  permutation  mappings  of  rectangular 
images  can  be  expressed  as  sequences  of  local  permutations,  or  transpositions,  in  the  image 
algebra  using  the  odd-even  transposition  sort  (GETS).  These  types  of  permutations  are 
required  for  implementing  the  DFT  locally. 

Definition  4.5.12.  A  permutation  template  t  is  an  template  with  the  properties  that  for 
every  x  €  X,  I71(x)l  =  1  and  if  y  €  T(x),  then  t^(y)  —  1. 

If  t  is  a  permutation  template  then  a^t  is  an  image  with  the  same  grey  values  as  a  but  (pos¬ 
sibly)  in  different  locations.  Henceforth  we  assume  that  X  is  an  m  x  n  array. 

Definition  4.5.13.  A  row  separated  permutation  template  is  a  permutation  template  t 


Up  = 


2 

n 

k— p— 1 


•w 


0,...,0,1  1  0 

0  0 


u  — 1  1 

.  0  -1 

.  .  0 


0  0  0 


0  . 

1  0 

-1  1 


*0 

Xq+Xi 

Xo— Xi 
X3-*2 

X4-X, 


*0 

Xo+Xi 
X0+X2 
X3-X2 
X.— X, 


*0 

Xo+Xi 

Xq-Hx, 

X0+X3 

*4-X3 


*0 

Xo+Xi 

Xq-HXo 

X0+X3 

X0+X4 


Figure  18.  Local  Algorithm  for  Computing  U5X 
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with  the  property  that  =  1  implies  that  u  =  x  . 

A  row  separated  template  will  permute  each  row  separately  and  simultaneously  but  possibly 
with  different  permutations.  These  types  of  permutations  are  actually  more  general  than  are 
required  for  the  local  DFT  but  the  GETS  algorithm  will  implement  them  just  as  easily  as 
the  more  restrictive  case.  If  t  is  as  in  Definition  4.5.13  then  the  mapping  t  >—*■  a0t  is  called 
a  row  separated  permutation.  The  same  definitions  and  remarks  can  be  made  for  the 
columns. 

The  GETS  applied  to  the  y*''*  row  of  an  image  is  given  by  the  following  rules: 
Alternate  the  following  steps: 

(1)  if  y  is  odd  and  a(x,y)  <  a(x,y+l),  then  switch  a(x,y)  and  a(x,y+l),  and 

(2)  if  y  is  even  and  a(x,y)  <  a(x,y+l),  then  switch  a(x,y)  and  a(x,y+l). 

The  output  will  be  an  image  with  each  row  sorted  in  increasing  order.  It  can  be  shown  that 
the  GETS  will  take  at  most  n  parallel  steps  to  execute  where  i.  and  ii.  each  count  as  one 
step,  Reference  25. 

We  now  write  the  GETS  in  the  image  algebra.  The  idea  is  that,  for  a  given  permu¬ 
tation  an  image  d  will  be  defined  containing  the  inverse  of  the  permutation  in  each  row.  This 
is  assumed  to  have  been  done.  The  sequence  of  transpositions  required  to  sort  d  will  be 
applied  to  an  input  image  a.  The  result  will  be  that  the  rows  of  a  are  permuted  by  the 
given  permutation. 

For  example,  if 


12  3  4 

3  12  0 

a  = 

5  6  7  8 

and  d  = 

0  2  3  1 

9  10  11  12 

0  2  13 

then  the  output  image  will  be 

4  2  3  1 

5  7  8  6 

.9  11  10  12 
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The  first  step  is  to  define  the  required  templates.  In  what  follows,  we  do  not  expli¬ 
citly  define  the  configurations.  If  a  template  in  question  has  nonzero  weights  at  a  point,  that 
is,  if  (x,y)  GX  such  that  there  is  a  (u,v)  €  X  with  ^x,y)(u>v)/0,  then  we  take  the  configuration 
at  (x,y)  to  be  all  such  points  (u,v).  If  not,  then  we  take  the  configuration  at  the  point  (x,y)  to 
be  the  set  {(x,y)}.  The  template  definitions  are  then  given  by: 
to: 


If  y  is  odd  and  y  ^  n-1,  then  for  every  x 


1  if  (u,v)  =  (x,y-Hl) 
0  else 


If  y  is  even  and  y  ^  0,  then  for  every  x 


^x.y)(u,v) 


1  if  {u,v)  =  (x,y-l) 
0  else 


te- 


Otherwise  for  every  x  and  {u,v),  t(^y)(u,v)  =  0. 


If  y  is  odd,  then  for  every  x 


1  if  (u,v)  =*  (x,y-l) 
0  else 


If  y  is  even  and  y  ^  n-1,  then 


*(x.y)('‘-’^) 


|l  if  (u,v)  =  (x,y-t-l) 
\o  else 


If  y  is  even  and  y  =  n-1,  then  for  every  x  and  (u,v),  =  0. 


If  y  is  odd  and  y  /  0,n-l,  then  for  every  x 


<y)(u.'^) 


1  if  (u,v)  =  {x,y) 

-1  if  (u,v)  =  (x,y-Hl) 
0  else 
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Otherwise  •(^)(u,v)  =  0. 


Se¬ 
lf  y  is  even  and  y  ^  n-1,  then  for  every  x 


«(xj)(«.v) 


1  if  (u,v)  =  (x,y) 

-1  if  (u,v)  =  (x,y-M) 
0  else 


Otherwise  for  every  x  and  (u,v),  «(^yj(u,v)  =  0. 
If  y  ^  0,  then  for  every  x 


’■(x,y)(u,v)  =  \ 


1  if  (u,v)  =  (x,y) 

1  if  (u.v)  =(x,y-l) 
0  else 


Otherwise  for  every  x  and  (u,v),  r(“  yj{u,v)  =0. 
If  y  ^  0,  then  for  every  x 


1  if  (u,v)  =  (x,y) 

1  if  (u,v)  =  (x,y-l) 
0  else 


If  y  =  0,  then  for  every  x 

|l  if  (u,v)  =  (x,y) 
=  \o  else 

Otherwise  for  every  x  and  (u,v),  r(^y)(u,v)  =0. 


The  templates  s,,,  and  r,,  are  depicted  in  Figure  19.  The  others  are  similar. 
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0 


1 


if  y  odd  and  y  #n— 1 


1  ‘ 

L2_ 

if  y  even  and 
otherwise 


if  y  odd  and  y  ^‘n— 1 
otherwise 


ify5>^o 


ify  =0 


Figure  19.  Templates  Used  to  Implement  GETS 
The  GETS  algorithm  takes  the  form: 

Input  image  a 

DOS  i  =  l,nmax 
DO  5  j  =  1,2 

IF  j  =  1  THEN  V  =  o  ELSE  v  =  e 

=  X>o(<l®«v) 
e  =  b©  Ty 
c  =  X<i(e) 
b  =  d*c+e*(d©tj 
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d  =  a*c-l-e*(a©tY) 
a  =  d 
d  =  b 


5  CONTINUE 

Output  the  permuted  image  a 


Note  that  due  to  the  use  of  the  symbol  =  to  denote  equality  rather  that  assignment 
we  do  not  write  statements  of  the  form  a  =  a*c.  The  parameter  nmax  can  be  set  to  the 
maximum  number  of  steps  that  the  GETS  will  take.  For  an  arbitrary  permutation  nmax  = 
n.  However,  if  n  =  mk,  where  m,k  >  1,  then  it  can  be  shown  that  the  GETS  will  execute  the 
generalized  shuffle  permutations  in  at  most  (m-l)(k-l)-fl  steps,  Reference  27.  Since  each  run 
through  the  inner  loop  represents  one  step  one  can  set  nmax  =  [(m— l)(k— l)/2l 

e.  Local  Implementation  of  the  100  x  100  DFT 

In  this  section  an  algorithm  for  the  implementation  of  a  100  x  100  DFT  on  a  mesh 
connected  array  using  the  image  algebra  is  developed.  We  are  assuming  that  the  array  is  at 
least  100  X  100.  As  pointed  out  earlier,  the  problem  is  to  derive  a  tridiagonal/permutation 
decomposition  of  the  matrix  F  jqq.  We  do  so  and  then  interpret  the  result  as  a 
template/image  decomposition.  We  write  the  algorithm  in  image  algebra  pseudo-code  and 
we  count  the  overall  number  of  parallel  steps  required.  We  assume  that  the  permutations  are 
precomputed  and  executed  by  the  routine  given  in  the  previous  section. 

Using  the  decomposition  theorems  the  following  equations  can  be  derived; 

There  exist  permutation  matrices  Qj,  Q2  such  that 

^100  “ 


Since 

F,  =  (F20l2)D(2,2)(l2®F2)P(2,2) 


it  is  clear  that 

F,;8)l25  =  P(.5O.2)(I^(8F2)P(2,5O)(D(2.2)«)Io5)(I.,<^P{25,2))(l5o0F.2)(Io(^P(2,25))(P(2.2)(^L5). 
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All  the  matrices  on  the  right  hand  side  of  the  last  equation  are  either  tridiagonal  or  permuta¬ 
tion  matrices.  As  for  the  other  factor 

F25  =  P(5,5)(l5<8)F5)P(5,5)D(5,5)(Is<8F5)P(5,5). 

Using  the  prime  decomposition  theorem,  there  exists  permutation  matrices  Rj,  R2  such  that 

F5=U5RiF11)A5F;(%oU5‘ 

where  the  a  needed  to  define  A  can  be  taken  to  be  2,  3,  or  4.  Continuing 
F1^>  =  (F2«)Ioyi^(2,2)(^)(Io(8)Fo)(^>P(2,2)t^>. 

F  urthermore 

(F2(8l2)^^^  =  P(2,2)(‘'(l2{8)F.2)(i’P(2,2)(^' 

Hence,  since  ,  we  can  put  these  equations  together  to  write  down  a 

tridiagonal/permutation  decomposition  of  I40F25.  Thus  we  arrive  at  a  decomposition  of  Fioq. 
In  order  to  compute  the  DFT  locally  along  each  row,  the  transformation  corresponding  to 
the  matrix  Iioo^Fioo  implemented  locally.  This  can  be  done  by  implementing  the 

sequence  of  transformations  indicated  in  Figure  20.  Since  all  matrices  appearing  in  the  table 
are  tensored  on  the  left  by  Ijqq  the  notation  is  suppressed.  It  can  be  seen  that  the  total 
number  of  parallel  steps  is,  at  most,  1014  permutation  steps,  52  addition  steps  with  at  most 
one  addition  per  pixel,  and  18  multiplication  steps  with  at  most  one  complex  multiplication 
per  pixel.  The  column  transforms  are  treated  similarly. 
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Matrix 

Type 

#  of  Parallel  Steps 

Operations/Pixel/Step 

l.)(l4®P(5,5))Q2 

P 

100 

2.)  l2o<8)U| 

A 

4 

1 

3.)Ioo®P(2,2)‘ 

P 

1 

4.)  I2O0R2 

P 

4 

5.)  l20®(l2®F2/'> 

A 

I 

1 

6.)  I2O0D(2,2)*(1>  or 

M 

1 

1 

6.)  Ioo<8)D(2,2f) 

M 

1 

1 

7.)  Ioo®P(2,2f ) 

P 

1 

8.)  Ioo0(l2®Fo)(‘) 

A 

1 

1 

9.)  l2o0P(2,2)(i> 

P 

1 

10.)  I2O0A5 

M 

1 

1 

11. )  Repeat  .S-9 

12. )  loo^Ri 

P 

5 

13.)  I2O0U5 

A 

4 

1 

14.)  I40D(5,5) 

M 

1 

1 

15.)  I40P(5,5) 

P 

17 

16. )  Repeat  2-15 

17. )  I40P(5,5) 

P 

17 

18.)  P(2, 2)01.5 

P 

100 

19.)  I20P(2,25) 

P 

24 

20.)  I5o0»^2 

A 

1 

1 

21.)  I20P(25,2) 

P 

24 

22.)  D(2, 2)0125 

M 

1 

1 

23.)  P(2,50) 

P 

100 

24.)  I5O0F2 

A 

1 

1 

25.)  QiPf50,2) 

P 

100 

Figure  20.  Sequences  of  Local  Transformations  Required  to  Compute  the  Row 
transforms  of  a  lOOX  100  DFT.  P  stands  for  addition  and  M  for  multiplication. 
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We  now  illustrate  how  to  use  the  operation  preserving  mapping  'if  ‘  to  write  the 
algorithm  in  the  image  algebra.  We  first  define  the  required  templates  and  then  the  required 
images  (recall  that  a  diagonal  matruc  corresponds  to  an  image).  In  these  template  definitions, 
we  take  the  configuration  at  a  point  (x,y)  to  be  the  set  of  all  points  (u,v)  such  that 

The  template  t  corresponding  to  Iioo<8l(Io<8lF2/^^  ^  given  by: 

If  y  mod  3=0,  then  for  every  x 

|l  if  (u,v)  =  (x,y) 

~  |o  else 

If  y  mod  3  =  1,  then  for  every  x 

1  if  (u,v)  =  (x,y) 

‘(x,y)('^,v)  =  1  if  (u,v)  =  (x,y+l). 

0  else 

If  y  mod  3=2,  then  for  every  x 

-1  if  (u,v)  =  (x,y) 

*(x.y)(u,v)  =1  if  (u,v)  =  (x,y-l). 

0  else 


The  template  s  corresponding  to  Iioo'^l5o'^f'2  i® 

If  y  is  even,  then  for  every  x 

1  if  (u,v)  =  (x,y) 

^(x,y)(“-'')  =  1  if  (u.v)  =  (x.y+1). 

0  else 

If  y  is  odd,  then  for  every  x 

-1  if  (u,v)  =  (x,y) 

«(x,y)(u,v)  =1  if  (u.v)  =  (x,y-l). 

0  else 

V 

The  templates  corresponding  to  the  decomposition  ol  lioo^iso^f^s  given  by 
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If  y  mod  5=0,  then  for  every  x 


ti: 


|l  if  (u,v)  =  (x,y) 
[O  else 


If  y  mod  5  =  1,  then  for  every  x 


1  if  (u,v)  =  (x,y) 

|1  if  (u,v)  =  (x,y-l). 
0  else 


Otherwise  for  every  x 


1  if  (u,v)  =  (x,y) 

-1  if  (u,v)  =  (x,y-l). 
0  else 


ti,i=  2,3,4: 

If  y  mod  5  =  i,  then  for  every  x 


1  if  (n,v)  =  (x,y) 

1  if  (u,v)  =(x,y-l). 
0  else 


otherwise  for  every  x 


1  if  (u,v)  =  (x,y) 
10  else 


ti©t2@t3©t4  will  correspond  to  Iioo^Ioo^^s-  templates  t[,  to,  t^,  corresponding  to 
defined  in  an  analogous  fashion.  The  preceding  templates  will  compute  all 
the  additions  required.  Note  that  the  configurations  associated  with  the  above  templates 
have  at  most  two  elements  each. 

The  images  corresponding  to  the  diagonal  matrices  are  now  defined.  These  images 
will  be  involved  in  computing  all  the  multiplications. 
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The  image  c  corresponding  to  Iioo^^so^^s  ^  defined  first.  Recalling  the  definition 
of  A5  it  can  be  seen  that 


if  y  mod  5  =  k  5^  0 

c  =  {  (x,y,c(x,y))  :  c(x,y)  =  1  ^  } 


where  \  =  f4(i'^  *),  i  =  (— 1)‘^.  Defining  c  in  terms  of  the  image  algebra  requires  no  if-then 
statements: 

4 

c  =  Sx-j(k  mod  5)*Xj 


where  k  is  the  image  used  to  define  the  DFT  template  and  \)  =  1. 

The  images  bj,  bo,  b3  corresponding  to  Iioo0D(2,2)(^Io5,  Iioo<8l40D(5,5),  and 
Ioo®D(2,2)^^)  respectively  are  now  defined.  Note  that  =-i  =  — {— 1)^/*.  Hence  D(2,2)  = 
diag(l,l,l,-i)  and  D(2,2)®Io5  =  blockdiag[Ios,Ios,L5,— ilosj.  Thus 

~  X>7s(^^)*(~*)  +  X<7s(^^) 


Let  u)  =  Wjs-  Observe  that 


D(5,5)  =  diag(l,l,l,l,l,l,w,w2^a^^a;^l,a^,w^a;®,u/,l,a^,u^,u^,w^2  1  . 


Hence,  defining  the  images  r  and  q  by  r  =  k  mod  5  and  q  =  INT((k-r)/5),  where  INT(x) 
denotes  the  integer  part  of  x,  we  have 


^>2  =  (<^25*I)**(r*q)- 


Similarly, 

ba  =  >U(r)*(-i)  +  X^{r) 

The  image  corresponding  to  the  conjugate  is  defined  by  replacing  i  with  -i. 

This  completes  the  definitions  of  the  images  and  templates  used  to  compute  the  row 
transforms.  Those  required  to  compute  the  column  transforms  can  easily  be  defined  by  inter¬ 
changing  the  roles  of  x  and  y  and  j  and  k.  The  code  for  the  column  transforms  can  also  be 
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integrated  into  the  code  for  the  row  transforms  by  adding  a  few  if-then-else  statements.  We 
now  write  the  adgorithm  in  the  image  algebra.  The  notation  PERMUTE(a,b)  denotes  a 
call  to  the  permutation  subroutine  where  the  image  a  is  permuted  and  the  result  is  put  into 
b.  Notice  that  due  to  the  use  of  the  symbol  =  to  denote  equality  rather  than  assignment  we 
do  not  write  statements  of  the  form  a  =  a*b. 

Input  100  X  100  image  a 

PERMUTE(a,b) 
a  =  b@ti©t3^©U@tJ 
PERMUTE(a,b) 

r  =  k  mod  5 
q  =  INT((k-r)/5) 

1*1  =  +  X<7s(^) 

bo  =  (W2s*l)**(r*q) 

1*3  =  X-4(**)*(-*)  + 

bs  =  X-4(*’)*>  +  W**) 

DO  10  m  =  1,2 
DO  5  n  =  1,2 
PERMUTE(b,a) 
b  =  a@t 
IF  n  =  I  THEN 
a  =  b*b3 
ELSE 
a  =  b*b3 
ENDIF 

PERMUTE(a,b) 

IF  n  =  1  THEN 
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a  =  b*c 


b  =  a 
ENDIF 

5  CONTINUE 

PERMUTE{b,a) 
b  = 

a  =  b*b2 

PERMUTE(a,b) 

10  CONTINUE 

PERMUTE(b,a) 
b  =  a@s 

PERMUTE(a,b) 
b  =  a*bi 

PERMUTE(b,a) 
b  =  a©8 

PERMUTE(b,a) 

Output  a,  the  row  transforms  are  complete. 

Algorithms  for  implementing  the  two  dimensional  DFT  of  any  size  can  be  derived 
using  this  methodology.  If  the  array  size  is  not  as  large  as  the  image  size  but  the  image 
dimensions  factor  into  products  of  integers  smaller  than  the  array  size  then  judicious  manipu 
lation  of  the  matrix  identities  can  lead  to  algorithms  for  breaking  the  computation  into  a  sue 
cession  of  smaller  two-dimensional  DFT’s.  For  example  suppose  a  1000  x  1000  DFT  is 
required  on  a  100  x  100  mesh  connected  array.  Then  since 

1000x1000  =  1000^^^  1000  ~ 

~  |F'(io^^ioo)^(^®’^^)(^io*^f^ioo)f*(^®’^^)|^  j'{Fio*8>Iioo)D(10,100)(l2Qi^F|QQ)P(10,100)  = 
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~  [(^  10®IlOo)®(F  10®  |i(IlO®  f' lOo)®  (^10®  ^ ioo)j 


*  (P(1(X),10)®P(100,10))  —  Pi(1 10000®  10x10)^2^^3(^100®  ^100x100)^4- 


where  the  P;  are  permutation  matrices  and  D  is  a  diagonal  matrix.  Thus  the  1000  x  1000  can 
be  broken  down  into  100  100  x  100  DFT’s  and  10,000  10  x  10  DFT’s.  One  could  also  use  25 
and  16  as  the  factors  rather  than  10  and  100.  In  any  case,  DFT’s  of  size  larger  than  the 
array  size  can  be  handled  in  this  way. 

By  the  convolution  theorem  this  technique  can  be  used  to  compute  any  circulant 
transform  locally.  Thus  the  methods  presented  here  offer  two  alternative  techniques  for  com¬ 
puting  separable  circulant  transforms  locally. 

f.  Local  Inverse  Templates 

We  now  present  an  example  of  how  Theorem  4.4.11  can  be  used  to  develop  parallel 
algorithms.  Let  X  be  an  8  x  8  array  and  let  m  €  circulant  local  average  tem¬ 

plate,  i.e.  fny(x)  —  if  X  is  an  eight  neighbor  of  y  and  0  else  and  m  wraps  around  at  the 

boundary.  It  can  be  verified  by  direct  calculation  that  m  is  invertible  with  circulant  inverse 
=  n  shown  below.  Thus  (a©m)©m~^  «  a.  In  particular,  we  can  apply  the  local 
averaging  transform  (a©m)  and  then  recover  the  original  image  a  by  applying  the  template 
to  the  averaged  image,  or  we  can  first  do  the  transform  a©m~^  to  obtain  the 
camouflaged  image  shown  in  the  bottom  left-hand  corner  in  Figure  21,  and  then  apply  the 
local  averaging  filter  to  this  transform  in  order  to  recover  a. 
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Figure  21.  Inverse  Averaging  Example 

The  upper  left  image  in  Figure  21  is  the  source  image  a;  the  upper  right  image  is 
the  averaged  image,  a@in;  the  lower  left  image  a©m~^;  and  the  lower  right  image  is 
(a@m)@in“^=(a@m“^)@m^,  that  is,  the  average  of  the  inverse  average  of  image  a. 

As  shown  in  Figure  22,  the  template  n  is  a  global  template  in  the  sense  that  the 
new  value  at  a  point  depends  on  the  value  at  every  other  point. 
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Figure  22.  Inverse  of  the  Local  Averaging  Filter  on  an  8X8  Array 

However,  n  is  separable  with  n  =3@t.  The  templates  s  and  t  are  shown  in  Figure  23 
below. 


Figure  23.  Separable  Components  of  ihe  Inverse  .A.veraging  Filter  n 

We  show  how  t  can  be  decomixeed  into  local  templates,  the  decomposition  of  s 
being  similar.  Sin''e  t  is  separable  and  circulant.  Theorem  4.4.11  applies.  The  point  0  is  a 
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minimal  point  for  t.  The  polynomial  representative  of  t  at  0  is  given  by 
Pt(x,y,0)  =  1— 2x4-x^-Hc®— 2x^+x^+x®— 2x^. 

The  roots  of  Pt(x,y,0)  are: 

qi  =  0.749927035578 

92,013  =  0.678637719963  ±  0.47903181019i 
94,95  =-1.00449681937  ±  0.3532353691i 
98,97  =  0.200895581622  ±  0.901038638456i 

where  i  =  (—1)^^^.  Hence,  if  a  is  an  arbitrary  image  on  X  then  a@t  can  be  computed  locally 
as  -2  *  (((((((a@qi)©q2)©<l3)  ")®®l7)>  where  the  qj  are  of  the  form  indicated  in  Figure  17. 

6.  LATTICE  STRUCTURES  IN  THE  IMAGE  ALGEBRA 

The  next  sections  discuss  the  lattice  substructures  inherent  in  the  Image  Algebra.  Inves¬ 
tigation  of  properties  involving  Q),  V,  and  0  has  lead  to  the  identification  of  lattice  struc¬ 
tures  and  some  useful  applications.  Discussed  here  are  the  vector  lattice  of  images  and  the 
structures  of  templates  under  ©  and  0  .  The  relation  of  mathematical  morphology  to  the 
Image  Algebra  is  discussed.  Properties  involving  images  and  templates  are  also  presented. 

a.  The  Vector  Lattice  of  Images 

Let  X  be  a  coordinate  set  with  a  finite  number  of  elements.  X  =  {  x^,  ...,  x„}. 
Define  V' :  (R)^  —*  R"  by  a  v^,  where  v,^  =(a(Xi),...,a(Xn)). 

The  function  i/^  is  an  isomorphism  from  R^  to  the  real  vector  space  R".  It  is  a 
trivial  matter  to  show  that  0  is  one-to-one  and  onto;  to  show  that  0(a  -I-  b)  =  0(a)  -I-  0(b), 
we  see  that  for  c  =  a  -I-  b,  0(c)  =  (a(xi)  -I-  b(xi),....a(Xn)  -I-  b(Xn))  =  -I-  v^  =  0(a)  -I-  0(b), 

using  the  definition  of  vector  addition.  The  remaining  properties,  namely  0(Q'a)  =  a0(a), 
a  €  R,  and  0(a«b)  =  0(a)*0(b)  are  similarly  easy  to  prove. 

Definition.  A  partly  ordered  vector  space  V  is  a  partially  ordered,  additive  (abelian)  group 
in  which  x  >  0  in  V  and  X  >  0  in  R  imply  Xx  >  0  in  V. 

Definition.  A  vector  lattice  is  a  partly  ordered  vector  space  which  is  also  an  1-group. 

For  example,  for  any  cardinal  number  n,  R"  is  a  vector  lattice.  Thus,  we  have  that 
is  a  vector  lattice,  and  any  results  known  about  vector  lattice  are  directly  applicable  to 
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the  corresponding  structure  of  images.  A  few  properties  are  listed  here. 

a  +  (b  V  c)  =  (a  +  b)  V  (a  +  c ) ,  V  a,  b,  c  G  .  (10) 

la  +  bl  <lal  +|bl  ,  where  I  a|  =  aVO  —  aAO  >  0  unless  a  =  0,  V  a,  b  GR^.  (11) 

7.  THE  ALGEBRA  DETERMINED  BY  ® 

a.  Operations  Between  Templates  under  Q 

Most  results  presented  here  are  for  the  general  case  of  templates  from  Y  to  X, 
though  of  course  they  hold  for  the  specific  case  when  Y  =  X.  We  will  also  discuss  the  con¬ 
cept  of  equivalence  classes  of  templates  which  is  used  to  describe  the  structures  identified. 

Let  i  G  R,  and  define  R^j^  =  {t  G  Ry-*x  •’  ^(x)  >  i  if x  €  T(y) }  and  R^g  = 

{a  G  R^ ;  a  >  0}. 

We  begin  with  a  basic  identity. 

Lemma  4.7.1.  s  ©  t  =  — (s  ®  (-^i)),  for  s  G  R^x  ’  ^  ^  • 

Proof;  Both  templates  are  members  of  Ry^x  >  same  configuration.  As  for  the 

gray  level  values,  on  the  left  hand  side  of  the  equality  the  image  at  location  y  G  Y 
has  form  {(x,rJx)):  ry(x)  ==  A  {ty(w)s^(x)  ;  x  G  S(w)}}.  On  the  right  hand  side 

^  w€T(y) 

of  the  equality  the  image  at  location  y  G  Y  has  form 

{(x,ti^(x)):it^(x)  =  -  V  {-f^(w)ajx)  ;  x  G  5(w)}}.  Since  in  general  -{  V  m„} 

w€T(y)  a€A 

=  A  {— nia  }  fo'’  indexing  set  A  and  a  collection  {mj,  }a€yi,  we  have 

a€A 

»'y(*)  =  A  { ty(w)«^(x)  ;  X  G  S(w)  }  =  -  V  {  -ty(w)8^(x)  :  x  G  5(w)  }  =  «  (x) 

w€T(y)  w€T(y)  ' 

Q.E.D. 

Lemma  4.7.2.  s®t  <  s®t  for  s  G  R^^x  ’  ^  ^  ■ 

Proof:  We  have  to  show  that  for  r  =s®t  and  for  u  =s®t  that  ry(x)  <  Uy(x),  as  the 
configurations  are  equal.  We  have: 

r  (x)  =  A  { ty(w)s^(x)  :  x  G  5(w)  }  and  since  all  values  are  positive,  the  product 
is  also  positive,  and  thus  we  know  A  { ty(w)«^(x)  :  x  G  5(w)  }  < 

w€  T(y) 

V  { ty(w)«  Jx)  ;  X  G  5(w)  }  =  Uy(x). 
wer(y) 
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Q.E.D. 

Lemma  4.7.3.  {  Ry_»x>  ^  }  is  a  poset. 

Proof;  We  use  the  definition  s  <  t  <=£>  ${y)  C  T(y)  and  <  ^,  V  y  €  Y.  This  obviously 
puts  a  partial  order  on  Ry-*Xi  ^  see  that  V  r,  s,  t  €  Ry-»Xi 

i.  t  <  t 

ii.  t  <  8  and  s  <  t  ^  a  =  t 

iii.  r<8,  s<t=s>r<t 

Q.E.D. 

That  not  every  pair  of  templates  are  related  can  be  seen  by  the  example  presented 
in  the  next  figure,  where  a,  t  €  Rx-*x  • 


1 

-3 

/  \ 

7 - ^ 

2 

3 

2 

a  = 

1 

\  / 

-1 

0 

Figure  24.  Two  Templates  Not  Related  Under  < 


Here,  a  ^  t,  t  ^  a. 

Theorem  4.7.4.  (s®t)®  r  =  s®  (t®  r)  for  r  €  R^^^  ,t  €  R^^y  ^  • 

Proof:  Let  u  =  s®t  €  R^^x  ’  ^  ~  (3®t)®  r  =  u®  r  6  Rw-x  -  let  q  =  t®  r  €  R^^y 
p  =  s®(t®r)  =s®q  €  R^Lx  •  f'i*‘st  we  show  that  the  configurations  K  and  P 
are  equal,  then  we  show  that  the  gray  value  functions  k  and  p  are  equal. 


<(w)  =  P{yv): 

X(w)=  U  U{z)=  U  f  U  S(y)l, 

i€;e(w)  i€;j(w)  U6T(i)  j 


and;i[w)=  U  S(f )  = 

r6s2(w)  fg 


U  S(f). 
U  r(g) 


g€*(w) 


Let  m  €  /C(w).  Then3  a  €>?(w)  such  that  m  E  U{z).  This  implies  thatB  a  €.*J(w) 
and  y  €  T(z)  such  that  m  E  5(y).  We  need  to  show  that3  f  E2(w)  such  that  m 
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G  S(f),  that  is,  that3  g  G,C(w)  and  f  GT(g)  such  that  m  G  5(f).  So  let 
g  =  *.  f  =y-  Then  g  G  £(w)  and  f  G  T{g)  where  m  G  S(f).  Thus,  X'(w)  C 
Conversely,  let  m  G  ^w).  Then 3  f  GQ(w)  such  that  m  G  5(f),  that  is, 

3  g  G  JZ(w}  and  f  G  T(g)  such  that  m  G  5(f).  Let  z  =  g,  and  y  =  f .  Then 
z  G  .^('w)  and  y  G  T(z)  where  m  G  5(y).  This  implies  that  m  G  JC(w).  Thus, 

.^w)  C  X’(w).  If  X’(w)  =  4>,  then  suppose  /(w)  #  (f).  Then  3  m  G  /RJw)  =s» 

3  g  G  X(w)  and  f  G  T(g)  such  that  m  G  5(f).  Letting  z  =  g  and  y  =  f  we  see  that 
m  G  K(w},  which  is  a  contradiction.  Tiius,  if  ^(w)  =  <^,  then  =  <^,  and 
KH  =  /(w). 

To  show  that  =  p„,  we  show  fc^(x)  <  Pw(x)  and  Pw(x)  <  fe*(x). 

fc^(x)  =  V  {  r^(z)  Uj(x)  :  z  G  ^(w)  and  x  G  U(z)  }.  Since  u^(x)  = 

V{iz(y)«y(x)  :  y  GT(z)  and  x  G  5(y)  },  we  have  fc*(x)  = 

V{»’,v(*)' [v{tz(y)«y(x);yGT(z)andxG5(y)}j:zG.?(w)  and  x€U(z)}  = 

V  {V  {r^(z)t^(y)ay(x)  :  y  G  T(z)  and  x  G  5(y)  }  :  z  G  ^(w)  and  x  G  U  5(v)}, 

ver(z) 

with  this  last  equality  following  from  the  fact  that  the  template  weights  are  positive. 
Let  D  = 

U  {  U  {r^(z)«z(y)«y(x)  :  y  G  T(z)  and  x  G  5(y)  }  :  z  G  ^(w)  and  x  G  U  5(v)}. 

v€  T(i) 

Next,  note  that  p^(x)  =  V  {  «i(x)  •  ^  €  2(w)  and  x  G  5(y) }  = 

V{  V{r.^^(n)t„(f):nG.?(vv)  and  fGT(n)}|s^x):fG  U  T(g)  and  x G  5(f) } 

J  ge^(w) 

=V{V{r.^(n)t„(f)«^x):nG.?(w)  and  fGT(n)};fG  U  7(g)  and  xG  5(f)}. 

gC^.fw) 

Let  E  = 

U  {U  {r.^^(n)t„(f)«^x)  :  n  G  /?(w)  and  fG  T(n)}  ;  f  G  U  7(g)  and  xG  5(f)}. 

g€ 

An  element  of  E  is  of  form  r.^,,(n)t„(f)«((x),  where  n  G  /?(w),  f  G  7(n),  and  x  G  5(f). 
We  will  now  show  that  D  C  E.  To  this  end.  choose  z  G  -?('w)  such  that  x  G  U(z}. 
Therefore, 3  y  G  7(z)  such  that  x  G  5(y).  Thus,  r.^^(z)tj(y)ay(x)  G  D.  To  show 
that  r.^^(z)tj(y)sy(x)  G  E,  we  must  show  that  for  n  =  z,3  f  G  ^(n)  such  that  f  =  y, 
and  that  x  G  5(f).  But  we  know  this  is  true,  by  choosing  f  =  y,  as 
y  G  7(z)  and  x  G  5(y).  Thus.  r^(z}t^{y)sy{x)  =  r.^^(nl<„{f)s,i(x)  G  E.  Therefore. 
A:^(x)  <  p^{x).  Conversely,  for  r^(n)t„(f)ar(x)  G  E.  where  f  G  7(n)  for  some 
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n  €  ^(w)  and  x  €  5(f),  we  must  show  that  r^(n)t„(f)«^x)  is  of  form  »‘^(a)ti(y)«y(x), 
for  some  z  €  ^(w),  where  y  €  T(z)  and  x  €  S(y).  Here,  we  choose 
z  =  n,  and  y  =  f,  and  we  are  done. 

Lastly,  we  note  that  if  K{-w)  =  then  Since  in  this  case  P[w)  =  <f),  we  have 

=  0  also. 

Q.E.D. 

Corollary  4.7.5.  (sOtlOr  =  8<2)(t(2)r)  for  r  €  ^  ^  • 

Proof:  The  proof  is  the  same  as  for  Theorem  4.7.4,  but  when  talcing  the  minimum  of  a  set, 
we  note  that  since  the  weights  of  the  templates  are  strictly  positive  on  their 
configurations,  the  minimum  is  never  zero. 

Note  that  if  one  of  s,  t,  or  r  is  not  positive,  then  Theorem  4.7.4  is  not  necessarily 
true.  For  example,  let  a,  t,  r  G  Rx-*X>  ^  Figure  25  below. 
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7 - ^ 

-3 

-9 

t  = 

1 

7 — ^ 
2 

r  = 

5 

7 - ^ 

25 

^ 

u=s<S)t  = 

-3 

7 - ^ 

-6 

-18 

qs=t©r  = 

5 

25 

7 - ^ 

50 

ii _ ^ 

^ 

k  =  u©r* 

-15 

-30 

'^-9o'1 

-450 

p  —8©  q  = 

-15 

-45 

7 - ^ 

-150 

-450 

^—2 

_ 1 

^ 

Figure  25.  (8©t)®r#8©(t®r) 


Thus,  obviously,  k  p. 

We  also  remark  that  since  the  template  weights  are  all  positive  in  Theorem  4.7.4 
and  Corollary  4.7.5,  that  an  abbreviated  definition  for  the  gray  value  function  of  r  =  8(St 
can  be  used  in  place  of  its  usual  definition  in  proving  these  results,  namely 

=  V{v(»)«*(x) :  *  e  7{y)}  , 

and  equivalently  for  ® .  This  is  because  if,  for  some  z  €  T{y)  there  is  an  x  ^  5(z),  then  •^(x) 
=  0,  which  makes  the  product  ty(z)«j(x)  =  0.  Since  all  such  products  are  positive  anyway, 
the  computation  of  rj,(x)  is  the  same  using  either  definition. 

Recalling  that  a  semigroup  is  a  set  having  an  associative  operation,  we  can  say  that 
Corollary  4.7.6.  {  ®  }  is  a  semi-group. 

Corollary  4.7.7.  {  ®  }  is  a  semi-group. 

If  {  Ry^x  I  V,  a  }  were  a  lattice,  we  could  obtain  a  richer,  more  useful  structure,  as 
many  properties  of  lattices  are  well-known.  However,  only  three  of  the  four  properties  neces¬ 
sary  (see  Reference  28,  for  example)  hold.  The  proofs  of  these  are  trivial: 

i.  8As=8,  8Va=8 

ii.  8At=tA8, 8Vt=tVs 
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iii.  (s  A  t)  A  u  *  8  A  (t  A  u),  (s  V  t)  V  u  =  8  V  (t  V  u), 

It  is  the  fourth  one  that  does  not  hold, 

iv.  8  A  (s  V  t)  =*  8  V  (s  A  t)  =  8 

since  by  definition  of  s  V  t,  the  resulting  configuration  at  y  is  S(y)  U  T(y),  which  is  not 
necessarily  the  same  as  S(y).  A  counterexample  is  easily  constructed.  This  is  presented  in 
Figure  26,  below. 


Figure  26.  s  A  (s  V  t)  ^  s 

That  8  A  (s  V  t)  8  follows  from  the  definition  of  equality  of  templates. 

This  makes  the  structures  (Ry^x  ’  ^  (^yIx  >  ^  ^  semi-lattices,  not  lattices. 

To  obtain  a  richer  structure,  we  use  the  equivalence  relation  between  templates  as  defined  by 
P.  Gader  6,  that  is, 

8  ~  t  if  and  only  if  ^  =  «y  V  y  €  Y. 

This  equivalence  relation  is  well  defined  under  ®  for  templates  and  images.  Thus  for 
Sj  ~  So ,  tj  ~  to ,  when  all  templates  are  elements  of  Rylx  ’ 

i.  [a,  ®  tj]  =  [so  ©  to]  =  fsil  ©  [til 
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ii.  [si  V  tj]  —  [s2  V  to]  —  [sil  V  [tj] 

iii.  a  Q  8]  =  a  ®  So,  if  a  >  0,  and  Sj  >  0. 

We  denote  this  set  of  equivalence  classes  by  E(R^^j. ).  Under  this  equivalence  relation, 
equality  of  two  templates  reduces  to  equality  of  the  two  gray  value  functions  only. 

Since  ~  is  well-defined  under  © ,  associativity  is  preserved  for  equivalence  classes  of 
templates,  that  is,  ([sj  ©  [t])  ©  [rj  =  [sj  ©  ([t|  ©  [r|).  Thus  we  have 

Theorem  4.7.8.  {  ),  ©  }  is  a  semi-group. 

Corollary  4.7.9.  {  ),  ®  }  is  a  semi-group. 

Lemma  4.7.10.  Let  [a],  [t|  6  EfR^^j^ ).  Then  (sj  A  ([sj  V  [tj)  =  [sj  V  ([a]  A  (tj)  =  [sj. 

Proof:  Since  [a]  =  {  u  €  Rylx  =  Ryl^  •  «y  =  V  y  €  Y  },  we  need 

only  to  show  that  the  gray  value  functions  are  equal,  as  mentioned  above.  Let  r  = 
a  V  t,  and  u=aAt,  p=aAr,  and  q  =  a  V  u.  Then 

Py(x)  =  «y(x)  A  {  ay(x)  V  ty(x)  },  and  qy(x)  =  fly(x)  V  {  Sy(x)  A  ty(x)  }.  Since  R 
itself  is  a  lattice,  we  see  immediately  that  Py(x)  =  qy(x)  =  «y(x). 

Q.E.D. 

This  helps  to  prove  that 
Theorem  4.7.11.  {  E(RYf,x ),  V  }  is  a  lattice. 

Proof:  This  follows  from  properties  i,  ii,  iii,  the  equivalence  relation  under 
©,  and  Lemma  4.7.10. 

Corollary  4.7.12.  {  E(Ry1x  ),  A  }  is  a  lattice. 

Much  is  known  about  lattices.  Some  properties  immediately  carry  over  to 
{  E(Ry^x  ),  V  }.  For  example,  we  can  say  that  {E(Ry^^  ),  V  }  is  a  distributive  lattice,  since 
R  is  itself  a  distributive  lattice.  This  follows  trivially  from  the  fact  that  the  gray  value  func¬ 
tion  of  s  V  t  is  computed  pointwise,  thus  inducing  the  distributive  lattice  property  of  the  real 
numbers  on  E(Ry^j^  )■  Also,  any  distributive  lattice  is  modular.  Thus,  we  can  apply  the  full 
theory  of  distributive  and  modular  lattices  to  E(Ry^y  )•  As  an  example,  in  any  distributive 
lattice  [sj  A  ([tj  V  (rj)  =([8]  A  [tj)V([s|  A  [rj).  There  are  too  many  theorems  to  state  them 
all  here,  but  Reference  28  can  be  consulted  for  more  results. 

A  multiplicative  semilattice  or  an  m-semilattice  is  a  semilattice  M  under  V  with  a 
multiplication  such  that  V  a,b,c  €  M, 
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a(b  V  c)  =  ab  V  ax:  and  (a  V  b)c  =  ac  V  c 


If  M  is  a  lattice  with  a  multiplication  and  the  above  property  is  true  also,  then  M  is  called  an 
m-lattice  or  an  l-groupoid.  An  m-lattice  which  is  a  semigroup  (monoid)  under  multiplication 
is  called  an  l-semigroup  (respectively  an  1-monoid,  short  for  lattice-ordered  monoid),  Refer¬ 
ence  28. 

A  multiplicative  poset  or  an  m-poset  is  a  poset  M  with  a  binary  multiplication 
which  satisfies 

a  <  b  implies  xa  <  xb  and  ax  <  bx 


for  all  a,  b,  X  G  M.  It  follows  easily  (using  the  above  property  of  an  m-semilattice)  that  any 
m-semilattice  is  an  m-poset.  We  will  show  that  some  of  these  substructures  exist  in  the 
Image  .\lgebra. 

Theorem  4.7.13.  s  ®  (t  V  r)  =  (s  ©  t)  V  (s  ©  r)  if  t,  r  G  ®  ^  ■ 

Proof:  Let  u=tVr,  k=8©  (tVr)=  3©u,  p  =  s©t,  q  =  s©r,  and  1  = 

(s  ©  t)  V  (s  ©  r)  =  p  V  q. 


ThenK(y)=  U  5(z)  =  U  S(a)  U 
*€T(y)U;e(y)  . 


u 

ni6^(y) 


5(m). 


Since  £(y)  =  ^y)U2(y)  =  U  5(z)U  U  5(m),  we  see  immediately  that 

z€T(y)  m€«(y) 

K(y}  =  L(y).  As  for  the  gray  values  fcj,(x)  and  ly{x}: 
ky(x)  =  V  {uy(z)s^(x)  :  z  G  U{y)  and  x  G  S{z)  } 

=  V  { [ty(*)  V  ry(z)l«j(x)  :  z  G  T(y)U^(y)  and  x  G  S{z)  } 

=  V  { ty(z)«i(x)  V  ry(z)«j(x)  :  z  G  T(y)U-?(y)  and  x  G  S(z)  }, 

with  this  last  equality  following  from  the  fact  that  the  template  weights  are  positive. 
Continuing, 

fcy(x)  =  V  {  fy(z)«j(x)  :  z  G  7(y )  and  x  G  S(z) }  VV  { ry(z)aj(x)  ;  z  G  /?(y)  and  x  G  5(z)  }, 
as  r,  s,  and  t  are  positive  and  have  zero  weight  value  outside  of  their  respective 
configurations.  This  is  exactly  ly{x).  Q.E.D. 

Corollary  4.7.14  (t  V  r)  ©  s  =  (t  ©  s)  V  (r  ©  s),  if  t,  r  G  ®  ^  ■ 

Proof:  We  omit  this,  as  it  is  very  similar  to  the  proof  of  Theorem  4.7.13. 
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Corollary  4.7.15.  a  ®  (t  A  r)  =(«  ®  t)  A(8  ®  r)  if  t,  r  €  ,  and  s  6  ■ 

Corollary  4.7.1ft  (t  A  r)  ®  a  =  (t  ®  a)  A  (r  ®  s),  if  t,  r  €  R^^  .  and  a  €  R^l^  . 

At  this  point,  a  number  of  results  follow  immediately  as  a  direct  consequence  of  the 
previous  properties.  For  example,  as  a  consequence  of  Theorems  4.7.8,  4.7.11,  4.7.13,  and 
Corollary  4.7.14,  we  have 

Theorem  4.7.17.  {  E(R^^^ ),  V,  ®  }  is  an  m-lattice. 
and  the  dual  result: 

Corollary  4.7.18.  {  E(Rx^3j ),  A,  ®  }  is  an  m-lattice. 

Lemma  4.7.19.  Let  [a),[t]  G  E(R2^y  )  such  that  [sj  <  [t],  and  let  [r]  G  E(R^^2  )> 

[ulG  E(R>°  ^ ).  Then  [u]  ©  [a]  ©  [r]  <  [u]  ©  [t]  ©  [rj. 

Corollary  4.7.20.  Let  [8l,[t|  G  E(R^®y  )  such  that  (si  <  [t],  and  let  [rj  G  ), 

[u]G  E(R>1^).  Then  [uj  ®  [a|  ®  [r]  <  [uj  ®  [t]  ®  [r]. 

Corollary  4.7.21.  ([a!  A  (tl)  ©  ([si  V  [tj)  <  ([t]  ©  [s])  V  ([s]  ©  [t])  for  [s],  [t]  G  E(R>°  ^ ). 
Theorem  4.7.22.  {  E(R^^^ ),  V,  ©  }  is  a  lattice-ordered  semi-group  (an  l-semigroup). 
Corollary  4.7.23.  {  E(R^f,jj ),  A,  ®  }  is  a  lattice-ordered  semi-group  (an  l-semigroup). 

Recalling  that  a  monoid  is  a  set  M  wl^ich  has  an  associative  binary  operation  and 
an  identity  element  I  G  M  such  that  I*x  =  x*I  =  x  V  x  G  M,  it  is  trivial  to  show  that 

Lemma  4.7.24.  {E(R^^jj )(  ©,  [I] }  is  a  monoid  where 

[Ij  =  {t  G  R^lx  •  *x(x)  *  1.  and  t^{y)  =  0  if  y  x  }. 

Lemma  4.7.25.  {E(R^^jj ),  ® ,  [I) }  is  a  monoid  where  I  is  as  in  Lemma  4.7.24. 

Theorem  4.7.17  and  Lemma  4.7.24  prove  that 
Theorem  4.7.26.  {  E(R^^x  ),  © ,  V,  [Ij  }  is  an  l-monoid, 
and  of  course  its  dual; 

Corollary  4.7.27.  {  E(R^^^ ),  ® ,  A,  [Ij  }  is  an  l-monoid. 

We  now  view  the  set  X  of  image  coordinates  as  a  rectangular  parallelepiped.  More 
specifically,  let  (nj,  nj,  nj^}  be  a  set  of  positive  integers  and  X  = 

{(x,,  ....  x^)  :  0  <  Xj  <  nj,  X;  G  Z^  i=l,  ...,  k  }. 
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Definition.  If  i  =  (zi,  •  •  •  ,Z|j)  €  Z*',  then  the  congruence  class  of  z  mod  X,  denoted  by 
<»>,  is  defined  as<z>  =(zi  mod(nj+l),  mod(nij4-l)).  If  Y  C  Z*'  then  the  congruence 
class  of  Y  mod  X  is  defined  as  <Y>  =  {<y>  :  y  €  Y}.  We  observe  that  <z>  €  X  <Y>  C  X 
<Z^>  =  X  =<X>,  and  if  x  €  X  then  <x>  =  x. 

Lemma  4.7.28.  If  t  is  translation  invariant,  then  t'^{y)  =  ty(x). 

Proof:  t',(y)  =  t*{<2x  -  y  >)  =  t^+^^{<2x  -  y  +  z>)-  ty{x),  where  z  =  y  -  x. 

Q.E.D. 

Lemma  4.7.29.  (s©t)y  =  tyQ>a' ,  a  €  Rx— x  >  ^  ^  ^y-Ix  ®  translation  invariant. 

Proof:  (8®t)y  =  {(z,  ry(z))  :  ry(z)  =  V{  ty(x)s,f(z)  :  x  G  T(y)  and  z  G  S(x)}  }  while 
^®a'  =  {(z,c(z))  :  c(z)  =  V{  ty(w)s'j(w)  :  w  G  5'(z)}  }.  We  have 

c(z)  =  V{^{w)s'j(w)  ;  w  G  S'Cz)}  =  V{yw)s'j(w)  :  w  G  T(y)} 

since  ty(w)  =  0  if  w  Tl(y)  and  Sy,  ty>Q. 

Continuing, 

c(z)  =  \/{  ty(w)s^(z)  :  w  G  T(y) }  =  V{  ^(w)s^(z) :  w  G  T(y)  and  z  G  5(w) }  since 
«^(z)  =  0  if  z  S(w).  This  is  exactly  ry(z). 

Q.E.D. 

Corollary  4.7.30.  (8®t)y  =  ty®8' ,  a  G  Rx^x  ’  ^  ^  x  ®  translation  invariant. 

The  next  two  results,  Theorems  4.7.31  and  4.7.32,  have  been  proven  in  Refer¬ 
ence  29.  We  shall  use  them  in  proving  Lemma  4.7.33. 

Theorem  4.7.31.  If  t  is  circulant,  then  t  is  translation  invariant. 

Theorem  4.7.32.  If  S  and  T  are  circulant,  then 

U  S(y)=  U  T(z)for  allxGX. 

y€T(x)  z€S(x) 

Lemma  4.7.33.  s®t  =  t®  s,  for  s,  t  G  Rxf.x  circulant. 

Proof:  By  Theorem  4.7.32,  the  configurations  are  the  same.  Thus  all  we  need  to  show  is 

(8®t)y  ==(t®8)y.  We  know  that  (8®t)y  =  ty®s',  as  by  Theorem  4.7.31  t  is  translar 
tion  invariant,  and  thus  the  gray  value  at  location  x  G  X  of  tyQa'  is  of  form  c(x) 

=  V  {  ty(w)s'Jw)  :  w  G  5'(x) } 
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■1 


=  V  {  <  w(y)»w{*)  •  *  ^  }  by  Lemma  4.7.28 

~y{ **<rw+x-w>(<y'+*~^>)«w(*)  •  *  €  5(w) },  by  the  invariance  of  t 

=  V  { t',(<y+x-w>)«^(x)  :  X  €  ^w) } 

=  V{  *  x(p)«<y+x-p>(*)  :  *  €  %+x-p) }.  letting  p  =  y  +  x  -  w; 

=*  V  {  *  x(p)*<y+x-p-Kp-x)>(<*‘*'P~*>)  -  P  S  %) }  by  the  invariance  of  s 
=  V{*'x(p)*y(p) :  P  ^  5(y)} 

—  y{  *  x(p)^(p)  ■  P  €  T'(x) }  as  the  product  is  zero  unless  p  €  T'(x), 
which  is  the  gray  value  at  location  x  of  Sy®t'  =  (t®8)y. 

Q.E.D. 

Corollary  4.7.34.  8®t  =  t(98,  for  s,  t  €  R-x^x  circulant. 

Theorem  4.7.35.  (  ),  ® ,  V,  [Ij  )  is  a  commutative  l-monoid,  where  C(R^^j. )  is 

the  set  of  equivalence  classes  containing  circulant  templates. 

Corollary  4.7.36.  (  ),  ® ,  A,  [I]  )  is  a  commutative  l-monoid. 

Several  comments  are  in  order  here.  First,  unlike  the  linear  algebra  structure  of 
the  Image  Algebra  which  was  shown  to  be  an  isomorphism,  we  have  no  isomorphism  between 
the  lattice  structure  and  subalgebras  of  the  Image  Algebra.  It  is  a  subtle  but  distinct 
difference.  Second,  the  structures  in  this  section  do  not  involve  images,  that  is,  elements  of 
R^  However,  some  useful  properties  of  images  and  templates  with  respect  to  ®  have  been 
identified,  and  these  are  present  in  the  next  section. 


b.  Operations  Between  Images  and  Templates  under  ® 

This  section  presents  properties  involving  images  and  templates  under  ® ,  V,  and-l-. 
There  is  no  well-known  structure  to  which  these  identities  adhere. 

Lemma  4.7.37.  (a  ®  s)  ®  t  =  a  ®  (s  ®  t),  for  a  €  R^g ,  s  €  R^^x  ’  t  €  Ry^.^  • 

Proof:  Let  c  =  (a® a)®t,  r  =  8®t,  and  d  =  a® (s®t)  =  a® r.  Since  both  c  and  d  are 

images  on  Y,  we  need  only  to  show  that  for  an  arbitrary  y  €  Y,  c(y)  =  d(y).  Note 
that  cW  =  ^^^a(xK(x))ya)  =  Also.  d(y)  - 

V  a(p)r  (p)  =  v  [\/{a(p)«y(n)«„(p)  :  n  €  T(y)  and  p  €  5(n) 

p€Xly)  il 


P6  y  :(in) ' 
meny) 
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U  U{a(p)^(n)«„(p)  :  n  €  T(y)  and  p  G  S{n)|.  An  eiement  of  Q  is  of  form 
p€  U  S(iii)l-  •* 

«erty| ' 

a(x)^(x)«j(x),  where  z  G  T(y)  and  x  G  S{z).  To  show  that  this  element  is  also  in  W, 

we  have  to  show  that  for  n  =  z  G  T(y)  that  3  p  G  S{n)  such  that  p  =  x.  Note  in 

enumerating  the  elements  in  W,  that  once  p  G  U  S(m)  is  chosen,  then  all  m 

mer[y) 

G  T(y)  such  that  p  G  S(m}  can  be  listed  and  the  inner  union  in  W  is  over  all  such 
m,  that  is:  for  p  G  U  S(m},  let  M  =  {mj,  .  .  .  ,  m^}  be  the  set  of  all  m  G  T(y) 

m€T{y) 

such  that  p  G  5(ni).  Then  U{a(p)^(n)a„(p)  :  n  G  T(y)  and  p  G  5(n)}  = 

k 

U  a(p)ty(mi)«„i(p).  Thus,  since  a(x)yz)«Jx)  G  Q  where  z  G  T(y)  and  x  G  S(z), 

set  p  =  X.  Hence  3  ntij  such  that  noj  =  z.  as  z  G  T(y)  and  p  =  x  G  S(z)  =  5(mi). 
Thus,  a(x)^(z)Sj(x)  =  a(p)^(mj)a,ni{p)  G  W  ^  Q  C  W.  Conversely,  choose  an 
element  a(p)ty(n)a„(p)  from  W.  We  must  show  that  3  z  G  T(y)  and  x  G  5(z)  such 
that  z  =  n  and  x  =  p.  Since  a(p)ty(n)s„(p)  G  W,  we  know  that  n  G  T(y)  where 
p  G  5(ii).  Simply  set  z  =  n,  x  =  p  and  we  are  done. 

Q.E.D. 

Corollary  4.7.38.  (a  ®  s)  O  t  =  a  ®  (a  O  t),  for  a  G  ,  s  G  ^  ,  and  t  G  R^l^  . 
Lemma  4.7.39.  a  ®  (s  V  t)  =  (a  ®  s)  V  (a  ®  t),  for  a  €R^o ,  s,  t  G  Rylx  • 

Proof:  Let  r  =  s  V  t,  b  =  a®  (s  V  t)  =  a®  r,  and  let  c  =  (a®8)  V  (a®t).  Let  y  G  Y.  We 
must  show  b(y)  =  c(y). 

We  have  b(y)  =  V  a|x)r  (x)  =  V  a(x)[B  (x)  V  ty(x)l  = 

.xe5(y)Ur(y) 


a(x)[s(x)Vty(x)jlV  V  a(z)(sy(z)Vty(z)jlV  V  a(w)[Sy{w)  V^w)] 

^x€5(y)\T(y)  j  ^xeny)Wy)  j  [w€S(y)nT(y) 

=  I  V  {a(x)«  (x)}|v  I  V  {a(z)ty(z)}lv  |  V  {a(w)«  (w)  V  a(w)ty(w)}l 

J  j  |w6S(y)nT(y)  j 

Note  that  of  these  three  terms,  the  last. 

V  {a(w)«j,(w)  V  a(w)tj,(w)},  can  be  written  as 
w6S(y)nT(y) 
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I  V  {a(w)«y(w)}  ,V  1  V  {a(w)iL(w)}]„  and  since 
|w€S(y)nT(y)  |'»€S(y)nT(y)  j 

V  {a(w)«y(w)}<  V  {a(w)«  (w)}and  V  {a(w)t^(w)}< 
w€5(y)nT(y)  »-€%)  w€S(y)nr(y) 

V  {a(w)i,(w)},  this  third  and  last  term  is  redundant  and  can  be  omitted.  Hence 
w€7Xy) 

Q.E.D. 

Corollary  4.7.40.  a  ®  (a  A  t)  =  (a  ®  a)  A  (a  ®  t),  for  a  GR^q  >  ^  ^yIx  • 

Lemma  4.7.41.  Let  a  €  R^g ,  t  €  Ry^x  •  Then  a©t=a©  t=a®  t  <=>  t  is  a  one 
point  template. 

Proof:  At  y  €  Y,  a©  t  has  gray  value  S  a(x)ty(x),  a©t  has  gray  value  V  a(x)t,(x), 

x€T(y)  ^  xer(y) 

and  a®t  has  gray  value  A  a(x)ty(x).  If  t  is  a  1-point  template,  then  obviously  all 

x€T(y) 

three  above  images  have  the  same  value,  namely  a(x)ty(x).  Conversely,  if 
S  a(x)ty(x)  =  V  a(x)IL(x)  =  A  a(x)ty(x),  and  if  T(y)  contains  more  than  one 

x€71(y)  '  x€T(y)  ^  x€7T[y) 

point,  we  will  not  necessarily  have  equality.  For  example,  choosing  t  6  R^lx  ^  ^ 
as  in  Figure  27, 

t  = 

Xj  Xo 

Figure  27.  Example  for  Proof  of  Lemma  4.7.41 

gives  us  a(xi)  +  ^{x^)  for  the  first  image,  V{a(xi),  a(xo)}  for  the  second  image,  and 
A{a(xi),  a(x2)}  for  the  third  image.  All  three  are  obviously  not  true  in  the  general 
case  for  every  a.  Choosing  a(xj)  =  1,  aJ{xn)  =  2  will  give  us  1  +  2  ^  2  5>^  1. 

Q.E.D. 

Lemma  4.7.42.  a  >  (a  ®  t~^)  ©  t',  where  a  €  R^g ,  t  €  R^lx  ’  ^  ^  T(y).  The  tem¬ 

plates  and  t'  are  defined  in  the  following  way:  t~^  =  (T,  «)  where 
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I  l/*x(y)  if  X  e  71(y) 

«y(x)  =  I  Q  i!  X  ^  T{y)  ■  general,  the  reflection  t'  is  defined  as  follows:  Let 

t  =  (T,  t)  €  Rx-*X-  Then  t'  =  (T  ',t')  is  defined  to  be 


T'(y)={z€X:z=2y-x,xGr(y)} 


and 


f  ty(<2y-x>)  X  €  r(y) 

*V(*)  =  I  0  X  ^  T(y)  • 

This  amounts  to  flipping  the  template  about  its  target  point.  An  example  is  given 
in  Figure  28. 


Figure  28.  Template  t  and  its  Reflection  t' 

Proof:  Let  x  €  X  and  write  7(x)  =  {Yq,  ■  ■  ■  ,  y„}  where  x  =  yo-  If  w  €  X  and  Z*'  —►Z*' 
is  a  translation,  then  when  writing  T(w)  =  (zq.  •  •  •  ,  *0  }>  '-h®  convention  that 

^(yj  =  *11  i^hat  is,  the  ordering  of  the  Zj’.s  with  respect  to  w  is  the  same  as  the  ord¬ 
ering  of  the  yj’s  with  respect  to  x.  We  shall  use  the  following  notation:  if 
T(yi)  =  {zq,  .  .  ■  ,  *„  }  and  yj  €  T(x),  write  Zj  =  yjj  and  ty^(Zj)  =  ty^(yij).  Thus, 

ty^lYok)  =  ^(yk)>  k  =0,...,n,  for  example.  Let  b  =  a®t“^  Then  b(yj)  = 

0  1 

A  {afy;:) — - — -},  i=0 . n.  Thus,  the  gray  value  at  location  x  of  (a® is 

j-o  *ri(yij) 

Myi)  <'yo(yoi)  =  *>(yk)  *'yo(yok)  =  Myk)  t'x(yk)  for  some  k.  Also,  we  know  that 
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“1  1 

**(yk)  =  *(ykj)".  X  =  a(yki)  T‘7 — r  for  some  i.  Since  €  T(x),  we  know 

that  X  G  T’(yij),  which  means  3  m  €  {  0,1,... ,n  }  such  that  x  =yijm-  Thus,  we  know 
that  ^^(ykm)  =  ^^(x)  =  t'x(yk)  f>y  properties  of  t'.  Also,  note  that 


a(ykm)  =  Kx)  =  a(yo)-  Since  b(y^)  =  a(yy)- 


1 


^v(yki) 

have,  for  the  gray  value  of  (a<2>t~^)©t'  at  location  x; 

1  a/  ^  \  1 


<  a(ykm)- 


^^(y  km) 


Myk)-*  x{yk)  =  a(yki)- 


*yk(yki) 


*  x(yk)  < 


^i((ykm) 


•tx(yk)  =a(x)- 


=  b(y J.  we 


■*'x(yk) 


=  a(x)  — — r'*x(yk)  ~  ^  desired  to  show. 

J  vlVu ) 


x(yk) 

Q.E.D. 

Corollary  4.7.43.  a  <  (a  ®  t')  ®  t“ where  a  €  ,  t  G  R^l^  ,  a^nd  y  G  T(y). 

The  proof  of  this  is  similar  to  the  proof  of  Lemma  -4.7.42  and  can  be  found  in  Reference  29. 
Lemma  4.7.44.  (a  +  b)(2it  <  a®  t  +  b®t,  for  a  G  ,  and  t  G  Ryl^ 

Proof:  Let  c  =  (a  +  b)®t,  and  let  d  =  a®t  +  b®  t.  We  must  show  that  for  y  G  Y  that 
c(y)  *  d(y). 

We  have  c(y)  =  =  ^gV^^[a(x)ty(x)  +  b(x)ty(x)l 


-  -  “(y). 


Q.E.D. 

Lemma  4.7.45.  (a  +  b)®t  <  a®t  +  b®t,  for  a  G  R^^ ,  and  t  G  Ryl^  ■ 

Lemma  4.7.46.  a® (t  V  s)  +  a©(t  As)  —  a@(t  4- s),  for  a  G  R^g ,  and  s,t  G  Ry^^  ■ 
Proof;  Let  b  =  a©  (tVs)  4-  a©  (t  As)  and  c  =  a©  (t  4-  s).  Then 

b(y)=  E  a(x)(ty(x)  V  3yix)|  4-  E  a(x)|  ty(x)  A  ej,(x))  = 

.x€T(y)US(y)  .x6T(y)US(y) 

E  a(x)[  ty(x)  V  ey(x)  4-  ^(x)  A  ey(x)  |  =  E  a(x)[  ty(x)  4-  ay(x)  ]  =  c(y). 
x€T(y)US(y)  x€T(y)USIy) 

Here,  we  use  the  fact  that  as  a  lattice  R  has  the  property  sVt4-sAt=s4-t, 

V  s,  t  G  R. 

Q.E.D. 
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Lemma  4.7.47.  (a@s)  V(a@t)  <  a@(8Vt),  for  t  €  ^  ^ 

Proof:  At  location  y  €  Y,  (a©8)\/{a©t)  has  gray  value 

I  S  I  5 

[x€S(y)  '  J  (*eT(y) 

Note  that  we  can  write  this  quantity  as 

I  E  a(w)8y(w)|  V  (  S  a(z)ty(z) 

(%»€S(y)UT(y)  j  [i€S(y)UT(y) 

as  if  *  ^  ${y),  then  ay(z)  =  0,  and  we  have  a(x)8y(x)  =0.  A  similar  case  holds  for 
any  z  ^  T{y).  Also,  on  the  right  hand  side  of  the  equation,  the  gray  value  is 

E  a(u)[8y(u)Vty(u)]. 
u€S(y)Ur(y) 

Note  that 

a(u)8y(u)  <  a(u)  |Sylu)  Vty(u)) 
a(u)<y(u)  <  a(u)  ( 8y(u)  V  ty(u)  ] 

as  a(u)  and  ty(u)  >  0,  This  implies  that 

E  a(u)8y(u)<  E  a(u)l«y(u)V^(u)) 
u€S(y)Ur(y)  u€5(y)UT(y) 

and 

E  a(u)ty(u)  <  E  a(u)[ay(u)Vty(u)j 
u65(y)UT(y)  u€r(y)UT(y) 

Thus, 

(  E  a(w)«y(w)l  V  I  E  a(z)ty(z)l<  E  a(u)[8^(u)Vfy(u)| 
[w€S(y)Ur(y)  J  \i€S(y)Ur(y)  J  ii65(y)UT(y) 

Q.E.D 

Lemma  4.7.48.  a  <  t  ^  a@8  <  a@t,  for  s,  t  €  ^  • 

Proof:  At  location  y  €  Y,  a@8  has  gray  value  E  a(x)ay(x).  Since  s  <  t,  we  know 

xesiy) 

%)C  T(y),  w  inch  implies  that 

E  a(x)«Jx)  =  E  aixlsJx)  +  E  aiz)sv(*).  which  is  true  as  a(z)aJz)  =  0 

x€S(y)  x€S(y)^  ieT(y)\S(y)^  ^ 


V  €  x€71[y)\5(y).  Thus,  since  «y  <  ^  V  y  €  Y, 

E  S  .a(x)<y(x)  +  S  a(x)ty(z)  = 

x€S(y)  '  i6T(y)\S(y)  '  *€S(y)  '  z€r(y)\S(y) 

S  a(x)^(x)  =  Zi  a(x)ty(x),  which  is  the  gray  value  of  a©t  at  location 
x€S(y)UT(y)  *€T(y) 

y€Y. 


Q.E.D. 

Corollary  4.7.49.  s  <  t  =s»  a©  s  <  a©  t,  for  s,  t  6  Ry—x  ^  ^  • 

Corollary  4.7.50.  a  <  t  =s*  a®  s  <  a®  t,  for  a,  t  €  Rylx  ^  ^  ^">0  • 


The  proofs  of  these  last  two  corollaries  are  similar  to  that  of  Lemma  4.7.48,  and  can  be 
found  in  Reference  29. 


c.  Applications  and  Examples  of  ©  Properties 


A  few  applications  of  the  preceding  results  are  now  presented.  Suppose  a  template 
t  can  be  decomposed  into  templates  t;  such  that 

t  =  tj  ©  to  ©  •  •  ■  ©  t„ . 

Then  a  ©  t  =  a  ©  (tj  ©  ©  •  •  •  ©  t„)  = 

((...((a  ©  t,)  ©  to  ©  )...)  ©  tj 

Computing  a  ©  t  using  the  above  equation  results,  in  most  cases,  in  fewer  multiplications 
than  computing  a  ©  t  directly. 

Now  suppose  one  cannot  find  templates  t-,  such  that 

t  =  tj  ©  to  ©  •  •  '  ©  t„, 

but  it  is  known  that  t  =  r  V  s,  where  r  and  s  can  be  decomposed  as  follows: 

r  =  rj  ©  ro  ®  •  ■  •  ®  r^ 

8  =  Si  ©  So  ©  •  •  •  ®  Sn, 

Then 


a  ©  t  =  a  ©  (r  V  s)  =  (a  ©  r)  V  (a  ®  s) 
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=  (a  ©  (r^  ®  Fj  ®  •  •  •  ©  r^))  V  (a  ©  (s^  ©  83  ©  •  •  ■  ©  s^)) 

=  ((...((a  ©  Fj)  ©  F2  ©  )...)  ©  fJ  V  ((...((a  ©  81)  ©  80  ©  )...)  ©  Sm) 
which  is  also  more  efficiently  computed  in  many  cases  than  with  t  directly. 

Here  is  an  example  of  this.  Let  t  €  Rx^x  template  as  presented  in  Figure 

29  below. 
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Figure  29.  Template  t 

Then  t  can  be  written  as  t  =  f  V  s  where  r  and  s  are  as  in  Figure  30. 
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Figure  30.  Templates  r  and  a 


Now  r  and  a  can  be  written  as  r  =  Tj  ®  To,  a  =  a^  ®  ao,  where  Tj, 
Figure  31. 


Figure  31.  Templates  Fj,  Fo,  Sj,  and  So 


Thus,  given  a  binary  image  in  which  we  wished  to  locate  an  object  which  at  least  is 
the  size  of  our  template,  we  could  replace  the  algorithm  a  t  by  the  last  expression  in  the 
equation 

a  ©  t  =  a  ©  (r  V  a)  =  (a  ©  r)  V  (a  ©  a)  = 

((a  ©  Ti)  ©  r.,)  V  ((a  ©  a^)  ©  So) 

The  number  of  multiplications  per  pixel  using  the  template  t  is  30,  while  the 
number  of  multiplications  per  pixel  using  the  templates  r  and  a  is  18. 

The  Division  Algorithm.  A  Euclidean  domain  is  defined  to  be  an  integral  domain  0  with 
a  Euclidean  valuation  on  D.  A  Euclidean  valuation  on  an  integral  domain  is  a  function  v  : 
D\{0}  — ►  {0,  1,  2,  •  ■  •  }  such  that  the  following  conditions  are  satisfied: 

(1)  For  all  a,  b  €  D  with  b  51^  0,  there  exist  q  and  r  in  D  such  that  a  =  qb  -I-  r,  where  either  r 
=  0  or  v(r)  <  t;(b). 

(2)  For  all  a,  b  €  D,  where  neither  a  nor  b  is  0,  ((a.)  <  oiab). 

For  example,  the  integers  Z  are  a  Euclidean  domain,  as  given  a,  t  6  Z,  there  exist 
unique  integers  q  and  r  such  that 

a  =  qt  +  r  ,  where  r  <  t. 

As  another  example,  take  the  polynomiab  in  1  variable.  Given  polynomials  a(x), 
t(x)  with  coefficients  in  a  field  F,  there  exist  unique  polynomials  q(x)  and  r(x)  such  that 
a(x)  =  q(x)t(x)  +  r(x)  ,  where  degree(r(x))  <  degree(t(x)). 

While  the  structure  of  templates  and  images  under  ©  and  V  comes  nowhere  close 
to  a  euclidean  domain,  it  has  been  possible  to  construct  a  division  algorithm  within  the  struc¬ 
tures  ©  and  EJ . 

For  our  purposes  we  will  use  a  valuation  function  defined  on  a  lattice  L,  Refer¬ 
ence  28. 

Definition.  A  valuation  on  a  lattice  L  is  a  real  valued  function  (functional)  u[x]  on  L  which 
satisfies 

i;(a)  4-  u(b)  =  i)(aAb)  +  u(a\/b). 

Given  an  image  a  and  a  template  t  we  can  define  a  division  algorithm  for  a  with 
basis  t  with  respect  to  ®  and  V.  We  define  a  valuation  function  on  the  lattice  (R^g ,  V)  by 
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i;(a)  =  Va  —  Aa. 

Let  E(P)  denote  the  sublattice  of  )  which  includes  only  invariant  templates, 

that  is,  E(P)  =  {  t  €  E(R^^3(  ) ;  t  is  invariant  }.  (Note  that  we  are  abusing  notation  here, 
and  using  t  €[t]  to  represent  [t],  an  equivalence  class.)  We  can  define  an  induced  valuation 
function  on  E{P)  by 

v:E(P)-*R, 

where  t  is  the  second  coordinate  of  t  =:  (T,t),  the  unique  extension  of  t  to  R".  Since  is  an 
image,  this  is  well  defined. 

We  present  a  few  examples  of  this  in  the  following  figure. 

5  3  9  6  9  9  7 

3  3  9  7  3  9  8 

2  9  1  6  4  9  4 

^^  2  297965 

5  2  9  5  8  9  7 

3  1  5  6  8  4  9 

t;(a)  =  8 

*  = 


0  0  0  0  0  0  0 

0  0  0  0  1  0  0 

0  0  0  1  19  I  0 

**^0000  1  00 
0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

t;(g=19 

Figure  32.  E.xamples  of  Valuation  Function 
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The  Division  Algorithm  for  ® .  Given  am  image  a  G  and  an  invariant  template 
t  G  Rxlx  ^  there  exist  images  q  and  r  such  that 

a  =  (q  ©  tO  V  r 

In  other  words,  an  expression  similar  to  the  division  algorithm  (a  =  q*t  +  r)  holds  for 
a  and  t.  Here,  we  have  «;((rit®  t~^)©  t^  <  t^t),  where  v  is  the  valuation  function  above. 
Thus,  in  this  sense,  the  remainder  r  is  small  with  respect  to  a  valuation  function. 

It  has  been  shown.  Reference  29,  that  choosing  q  =  a  ®  t~^,  and  r  = 
a  *  [x>o(a~(a  ®  t~\)  ©  t*)]  will  make  the  above  equality  true. 

One  useful  result  of  being  able  to  express  a  in  this  way  is  in  the  following  applica¬ 
tion.  Making  the  following  recursive  definitions,  we  can  write  a  in  a  series  form 

a  =  ro  V  (ri  ©  to  V  (r,  ©  V  •  •  •  V  (r„_i  ©  (t')""')  V  (a„  ©  (t')") 

=  {Vr,®(t')''}v(a.ffl(tT) 


where 

ao=». 

rk  =  aj.  *  [x>o(»k“(ak  ®  O  ®  ^')]  ,k  =  0,...,n-l, 

ak  =^-i  ®  =  l,...,n, 

and 

(t')*' =  t' ©  t' •  •  ■  ©  t',  k  times. 

An  Image  Blurring  Technique.  During  our  investigation  concerning  the  properties  © , 
several  novel  blurring  techniques  which  emphasize  hot  targets  and  blurs  streaking  in  infrared 
images  were  discovered.  Here  we  present  one  of  these  blurring  transforms. 

The  Image  Algebra  equation  is 

[(a©ti— a®to)"  -|-(a©t3— a®t4)“l'‘ 
with  the  templates  tj’s  as  below. 
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A  pictorial  example  is  given  in  the  next  two  figures. 


Figure  33.  Input  Image  a 
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Figure  34.  Blurred  (Transformed)  Image 

8.  THE  ALGEBRA  DETERMINED  BY  0 

The  algebraic  structure  determined  by  0  is  mathematically  not  as  rich  as  the  one 
under  ® ,  but  it  does  yield  some  surprising  results,  such  as  a  division  algorithm  similar  to  the 
one  for  © . 

a.  Operations  Between  Templates  under  0 

As  was  for  the  case  of  the  algebra  determined  by  © ,  the  results  are  presented  for 
the  most  general  situation  possible,  that  is,  for  templates  from  Y  to  X 

Two  very  basic  results  are  Lemmas  4.8.1  and  4.8.2: 

Lemma  4.8.1.  s  Q  t  =  -{—8  0  — t)  for  s  €  R^^x  >  ^  ^  ®'y1z  • 

Proof:  Both  templates  are  members  of  Ry^x  ’  same  configuration.  We  note 

that  since  the  template  weights  are  added  and  not  multiplied,  as  in  the  analogous 
case  for  ©,  that  we  introduce  a  minus  sign  on  both  of  the  templates  on  the  right 
hand  side.  On  the  left  hand  side  of  the  equality  the  image  at  location  y  G  Y  has 
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{(x,ry(x)) ;  ry(x)  =  y  w)  4-  •  Jx) :  x  €  5(w)  }  } . 

On  the  right  hand  side  of  the  equality  the  ima^e  at  location  y  €  Y  has  form 
{(x,tty(x));  tty(x)  ^  ^w)})}- 

Similar  to  the  reasoning  in  the  proof  of  Lemma  4.7.1  we  can  conclude  that 

V*)  =  ^  +  -«Jx)I  :  X  €  S{w)  }  ~ 

w€T(y) 

-V  +  -«w(x)  :  X  €  S(w)}  =  ■u^(x). 

w€T(y) 


Lemma  4.8.2.  s|2|t<s|2tfors€  R^x  >  *  ^  ®’y~z  • 


Q.E.D. 


Proof:  We  have  to  show  that  for  r  =  s  Q  t  and  for  u  =  s  Q  t  that  ry(x)  <  tty(x).  We  have 
weT(y) 


w€T(y) 


Q.E.D. 


We  also  have  the  associativity  of  templates  under  Q  : 

Theorem  4.8.3.  (s  Q  t)  ta  r  =  s  0  (t  Q  r)  for  r  €  R^^2  ,t  G  Rz.1y  >*  ^  ®'y~x  • 

Proof:  The  proof  of  this  is  very  similar  to  that  of  Theorem  4.7.4,  and  we  omit  it  from  this 
report.  However,  it  is  included  in  Reference  29. 

Corollary  4.8.4.  (s  □  t)  Q  r  =  s  Q  (t  Q  r)  for  r  G  Rwl.z  ^  ^-^^y  >*  ^  ^y^x  • 

Thus  we  know 

Corollary  4.8.5.  {  Rxlx  >  M  }  is  semigroup. 

Corollary  4.8.8.  {  Rx.lx  ’  23  }  is  a  semigroup. 


Similar  to  the  structure  under  ® ,  we  have  the  distributivity  of  ®  over  V,  and  their 


duals: 


Lemma  4.8.7.  s  0  (t  V  r)  =  (s  0  t)  V  (s  0  r)  if  t.  r  G  RyIz  .  s  G  R^^x  • 

Proof:  Let  u=tVr,  k=8  0  (tVr)=80u,  p=80t,  q=80r,  and  1  = 

(a  0  t)  V  (a  0  r)  =  p  V  q.  The  proof  that  the  configurations  are  the  same  is 
identical  to  that  found  in  Theorem  4.7.13.  As  for  the  gray  values  A^(x)  and  ly,(x), 


we  have: 

ky{x)  =  V  {  tty(»)  +  «,(x) ;  2  €  2i(y)  and  x  €  5(z)  } 

=  V  { [^(*)  V  ry(*)]  +  e,(x)  :  z  €  T(y)U^(y)  and  x  €  S(z)  } 

=  V  {(^(z)  +  «,(x))  V  (ry(z)  +  «i(x)) :  z  G  T(y)U^(y)  and  x  €  S(z)  }, 

with  this  last  equality  following  from  the  fact  that  the  template  weights  are  positive. 
Continuing, 

*^(x)  =  V  {  ^(z)  +  «,(x)  :  z  €  T{y)  and  x  G  5(z) }  VV  {  ry(z)  +  Sj(x)  :  z  G  ^(y)  and  x  G  5(z)  }, 
as  T,  s,  and  t  axe  positive  and  have  zero  weight  value  outside  of  their  respective 
configurations.  This  is  exactly  ^(x). 

Q.E.D. 

Corollary  4.8.8.  (t  V  r)  M  a  =  (t  0  s)  V  (r  Q  s)  if  t,  r  G  ^  ■ 

Corollary  4.8.9.  a  ES  (t  A  r)  =  (a  E3  t)  A  (s  0  r)  if  t,  r  G  Rylz  >  and  a  G  Rz^x  • 

Corollary  4.8.10.  (t  A  r)  0  a  =  (t  0  a)  A  (r  0  s)  if  t,  r  G  ®  ^ 

By  Lemma  4.7.3  and  Corollaries  4.8.7  and  4.8.8,  we  have 
Theorem  4.8.11.  {  R^lx  >  V,  0  }  is  am  m-semilattice. 

Corollary  4.8.12.  {  R^lx  >  ^  A  0  }  is  an  m-semilattice. 

Since  any  m-semilattice  is  an  m-poset,  we  have 
Theorem  4.8.13.  {  Rx^x  >  ^  V,  0  }  is  an  m-poset. 

Corollauo^  4.8.14.  {  Rxlx  ■  ^  A  0  }  is  an  m-poset. 

Thus  we  know  that  the  following  two  lemmas  are  true,  being  properties  of  an  m- 

poset. 

Lemma  4.8.15.  u<t*^s0u<s0t  and  u  0  a  <  1 0  a  for  a,  t  G  E(R2^y  ),  and 
r  S  E(RJ1, ). 

Lemma  4.8.16.  u<t^a0u<a0t  and  u  0  a  <  1 0  a  for  a,  t  G  E(R2^y  )- 
r  G  E(R>lx  )■ 

Lemma  4.8.17  u  a0u0r<s0t0rfor 

rGR^!.,,tGR^Y.seR?lx 
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By  Corollary  4.8.5,  we  have 


Theorem  4.8.18.  {  S  V  •  0  }  is  a  po-semigroup  (partially  ordered  semigroup). 

Corollary  4.8.19.  {  ^  Q  }  is  a  po-semigroup. 

It  is  not  obvious  how  to  put  a  nice  equivalence  relation  on  R>°^  to  make  it  a  lat¬ 
tice.  In  the  case  of  ® ,  the  induced  operation  ®  between  equivalence  classes,  [sj®  [t]  = 
[8®t],  was  well-defined.  However,  using  this  same  equivalence  relation,  we  can  find  exam¬ 
ples  of  s  and  t  where  [sj  |2  [t]  ^  [a  Q  tj.  Given  the  templates  s,  u,  and  t  as  in  Figure  35,  and 
the  equivalence  relation  ~  as  defined  in  Section  7  we  see  that  aE3t=s|2uis  not  neces¬ 
sarily  true  for  all  templates,  where  t  and  u  are  in  the  same  equivalence  class.  In  Figure  35, 
we  present  an  example  of  this. 


Figure  35.  Example  of  [sj  Q  |tj  [s  Q  uj,  t~u 

Here  we  see  that  for  t  and  u  in  the  same  equivalence  class,  s  Q  t  7^  s  0  u.  Thus, 
using  the  equivalence  relation  ~  and  the  naturally  induced  operation  0  between  equivalence 
classes,  we  find  that  this  induced  operation  is  not  well  defined.  There  may  be  another 
equivalence  relation  that  may  induce  a  well-defined  0  operation  between  equivalence 
classes,  but  it  is  not  as  readily  apparent.  Not  having  the  strength  of  a  lattice  structure,  but 
only  of  a  semi-lattice,  the  structure  of  the  Image  Algebra  under  0  is  not  as  rich.  We  state 
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a  few  more  theorems  which  may  prove  to  be  useful  at  some  future  date.  Any  proofs  not 
given  can  be  found  in  Reference  29. 

Theorem  4.8.20.  a  Q  t  =  (t  Q  s)  for  a,t  €  R-x— x  ’  ^  circulaat. 

Corollary  4.8.21.  a  E]  t  =  (t  E3  a)  for  s.t  €  Rx— x  >  *»  ^  circulant. 

Theorem  4.8.22.  {  Rxlx  >  ^  Q  ,  1  }  is  a  po-monoid,  where  i(x)  =  {x}and  tx{x)  =  1. 
Corollary  4.8.23  {  Rxl.x  >  ^  A,  Q  ,  i  }  is  a  po-monoid,  where  i  is  as  in  Theorem  4.8.22. 

b.  Image-Template  Operations  Under  Q 

As  was  the  case  for  image-template  operations  under  ® ,  there  Ls  no  well-known 
structure  for  the  image-template  operations  under  EJ .  Some  basic  properties  that  have  been 
found  to  hold  are  listed  here. 


Lemma  4.8.24.  a  □  a  =  -{—a  Q  -a),  for  a  E  R^,  a  €  Ry— x 


Proof:  a  □  a  =  {(z,  c(z))  :  c(z)  =  A  a(y)  -I-  a  (y),  y  6  Y},  and  -(—a  Q  -a)  = 

y€S(z) 

{(z,d(z)  :  d{z)  =  V  -a(y)  +  -«i(y),  y  G  Y}.  We  have 

y6S(z) 

*^(»)  =  V  -  Wy)  +  «z(y)l  =  -  A  [a(y)  -f  «,(y)l. 

y€S(*)  y€S(t) 

This  means  that  — d(z)  =  ^[a(y)  -f-  «j(y)]  «  c(z). 


Q.E.D. 


Lemma  4.8.25.  (a  0  a)  ia  t  =  a  0  (a  Q  t)  for  a  €  R^,  t  €  Ry— z>  *  G  Rz_»x- 

Proof:  Let  b  be  the  image  b  =  (a  0  a)  0  t.  and  let  c  =  a  0  (a  0  t).  Then 

My)  =  V  {[  V  a(x) -t- «j(x)] -1- lL(z)l=  V  I  V  a(x) -I- a,(x)  +  L(z)l  = 
i€rfy)|  U€S(ir  'j  i6T(y)|.x65Ur 

V  {a(x)  4-  «j(x)  -f  ty(z)  :  x€5(z) }.  Also,  we  know  that 
i€T(y) 


c(y)  =  V  Ja(w)  -4-  [  V  {ty(p)  a  (w)  :  w65(p)} 
we  U  5(u)  lpeT(y) 

ueT(y)  ^ 

V  I  V  {a(w)  4- (^(p)  4- «p(w)  :  w€5fp)}L  This  last  equality  is  true 

w€  U  s(u)  P^TIy) 

u€T(y)  '■  ' 

because  if  w  ^  ^(p)  then  a(w)  is  not  added  to  ^(p)  4-  «p(w).  Let  B  = 

U  {a(x)  4-  «i(x)  4-  ty(*)  ■  xG5(z)},  and  let  C  = 

*€T(y) 
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.M  I  U  +  *y(p)  +  «p(w)  :  we5(p) }  I  Choose  an  element  in  B.  It  is 

w€  U  s(u)  P^ny)  '  ( 

of  form  a(x)  +  «j(x)  +  Uz),  where  z  €  T(y)  and  x  €  S(z}.  Since  S(z)  C  U  S(u), 

u€r(y) 

we  have  x  €  U  5(u).  Let  w  =  x  and  p  =  z.  Then  a(x)  4-  sJx)  +  iLfz)  = 

u6Tl(y)  ' 

a(w)  4-  «p(w)  4*  ty(p)  €  C.  Therefore,  B  C  C.  Conversely,  choose  an  element 
a(w)  4-  ty(p)  4-  in  C.  Then  p  €T{y)  and  w  €  S(p).  Let  x  =  w  and  z  =  p. 
Then  a(w)  4-  ^(p)  4-  «p('w)  =  a(x)  4-  s^(x)  4-  ^(z)  G  B  and  thus  C  C  B.  Therefore, 
C  =  B,  and  we  are  done. 

Q.E.D. 

Corollary  4.8.26.  (a  0  a)  Q  t  =  a  E3  (a  0  t),  for  a  6  R^,  t  G  Ry~*Z’  ®  ^  ^-z-^x- 

Theorem  4.8.27.  (a  0  —a)  0  —fc  =  a  0  [—{a  Q  t)]  for  a  G  R^,  t  €  Ry-.z>  ®  ^  ^■z-.x- 

Proof:  By  Corollary  4.8.26,  we  know  (a0  —a)  0  — t  =  a0  (—a  0  — t),  and  by  Lemma 
4.8.1,  we  know  — a  0  — t  =  —(a  Q  t).  Thus,  (a  0  —a)  0  — t  =  a  0  [-{a  M  t)]. 

Q.E.D. 

Theorem  4.8.28.  (a  0  a)  Q  t  <  (a  0 1)  0  a,  a,  t  circulant,  and  a  €  R^,  t  G  Ry-.z» 
s  €  ®'Z-*X- 

Proof:  Let  d  =  (a  0  a)  0 1  and  let  e  =  (a  0 1)  0  a.  The  gray  values  of  d  and  e  at  y  G  Y 
au’e,  respectively, 

«i(y)  =  V  I  A  {a(x)  +  «Jx)  +  ty(z)} 
ieT(y)U€S(i) 

«(y)  =  J),  i  V  {a(u)  4-  Uu)  4-  «^(w)} 

w€5(y)ju€T(w)  ' 

Let  T(y)  =  {zj,  .  .  .  ,z„ },  and  let  S(zi)  =  {x-.i,  .  .  .  ,Xjk }.  Also  let 
S(y)  =  {wi,  .  .  .  ,w^  },  and  let  r(wj)  =  {u^,  .  .  .  ,Ui„  }.  Note  that  I  5(p)l  =  k,  and 
inp)!""  V  p  G  X,  as  a  and  t  are  circulant  templates.  Now  we  can  write  the 
gray  values  of  d  and  e  as 

d(y)  =  **i(*ii) 
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«(y)  =  {qi  W«p<,)  +  *w/Upq)  +  «y(Wp)  }| 

Letting  =  a(Xij)  +  s^^Xy)  +  ty(Zi)  and  fp,  =  a(Up„)  +  t^^(upq)  +  «y(Wp),  we  can 
write 

Since  I  {5,:}  I  l{^Da}l  ,  (as  the  templates  are  circulant),  if  we  can  show 

*  1»1  J*l  P““l  ' 

that  {&,j}  C  {Cpq},  then  we  will  have  {&,j}  =  {Cp^}.  Using  a  theorem  of  BirkhofTs 

called  the  min-max  property,  which  states  \)  Xa,i  <  X  'Da;:,  Reference  28,  we  will 

j-U-l  ^  i-lj-l  ^ 

have  our  result  as  desired.  To  this  end,  let  =  a(xij)  +  «i,(Xjj)  +  ^(Zj).  Since  s  is 
circulant,  it  is  translation  invariant,  and  we  have  y  =  <Zi  -H  v>,  for  some  v  G  Z^. 
Thus  <y  —  v>  =  Zj.  Let  m  =  <Xij  +  v>.  Note  that  <y  +  Xjj  +  v  —  y>  =  <Xij  +  v> 
=  m.  Therefore,  ^  -  y >(<*i  +  Xy  +  v  -  y>)  = 

t J<Zi  +  Xij  +  V  -  y>)  =  t„(<y  -  V  +  xy  +  V  -  y>)  =  t J<Xij>)  =  t Jxy).  Since  t 
is  circulant,  we  have  Xy  €  T(m).  Also,  e^^Cxy)  =  «<i,  +  v>(<Xy  +  v>)  =  «y(m). 

Again,  s  circulant  gives  us  that  m  €  5{y).  Therefore, 

a(Xij)  +  «i,(Xij)  +  *y(*i)  =  a(xy)  +  ijxy)  +  «y(in),  where  m  G  ${y)  and  Xy  G  T(m). 
Letting  Wp  =  m,  we  have  for  some  q,  Up^  =  Xy.  Thus 

^Xy)  +  *in(Xij)  +  «y(ni)  =  a(upq)  +  t^^(Upq)  +  «y(Wp).  Thus  we  know  {(Sy}  C  {fp^}. 
Applying  the  min-max  property,  we  are  done. 

Q.E.D. 

Lemma  4.8.29.  (a  Q  s)  Q  t  =  (a  Q  t)  (3  s,  s,  t  circulant  and  a  G  R^, 
s,t  G  Rx-*x- 

Proof:  We  know  that  (aI3  s)  0 1  =  al3  (s  E3  t).  Since  s  and  t  are  circulant,  we  have 
s  13  t  =  t  Q  s.  Thus,  a  13  (s  13  t)  =  a  13  (t  13  s). 

Q.E.D. 
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Corollary  4.8.30.  (a  Q  s)  t  =  (a  Q  t)  Q  s,  s,  t  circulant  and  a  6  R^, 
s,t  G  Rx_»x- 

Lemma  4.8.31.  (a  V  b)  Q  t  =  (a  M  t)  V  (b  0 1)  a  G  R^,  t  G  Ry-*x- 

Proof;  On  the  left  hand  side  of  the  equation  we  have  for  the  gray  value  at  y  G  Y, 

V  [a(x)\/b(x)j  +  ty(x)  , 
x€T(y) 

and  on  the  right  hand  side  we  have 


But  V  [a(x)Vb(x)l  +  ty(x)  =  J  V  (a(x)  +  tj,(x)llv  J  V  [b(x)  +  ly(x))i  = 


x6T(y) 


xeT(y) 


x6T(y)‘ 


+ ‘.W}l  V  +  y,)i). 


Q.E.D. 


Corollary  4.8.32.  (a  A  b)  Q  t  =  (a  Q  t)  A  (b  Q  t)  a  G  R^,  t  G  Ry-.x- 

Lemma  4.8.33.  a  >  (a  Q  — t)  Q  t',  where  a  G  R^q  ,  t  G  R^lx  ’  ^  reflection  of  t. 

The  proof  of  this  is  similar  to  the  case  involving  ©  (Lemma  4.7.42). 

Corollary  4.8.34.  a  <  (a  0  t')  |S|  —t, where  a  G  R^g ,  t  G  R|1x  • 

These  next  few  properties  involve  -H,  V,  and  0  . 


Lemma  4.8.35.  -^(a  0  s  +  a  0  t)  <  a  0  (s  + 1),  a  G  R^,  s,  t  G  R^l,, 


Proof;  At  location  yGY,  a08+a0t  has  gray  value 


and  a  0  (s  + 1)  has  gray  value 

V  a(w)  +  [  «j.(w)  4-  tj.(w)  1 
wer(y)USfy) 

Now  V  a(w)  4-  { Sj.{w)  4-  ty(w)  j  >  V  a(w)  4-  [  «y(w)  j  > 
w€T(y)US(y)  w6T(y)US(y) 

V  a(w)  4-  [«y(w)|,  as  «y(w)  >0  and  5(y)  C  T(y)US(y).  Similarlv, 

w5(y) 

V  a(w)  4-  ^  V  a(w)  4-  [  tv('''^)!-  Thus, 

w6T(y)U5(y)  ’ 
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where 

=  a, 

fk  =ak  *  !X>o(ak-[(ak  0  -t)  13  t'])),k  =0,...,n-l. 

^k  =  ^k-i  Q  -t.k  = 
and 
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(t')*'  St'  n  t'  -  •  •  H  t',  k  times. 

Note  that  must  be  non-negative  for  each  k=l,...,n,  in  order  to  continue  to  apply 
the  division  algorithm  to  it. 

It  can  be  shown,  Reference  29,  that  for  each  k  =0,1,2,..., 
v((rk  E3  -t)  0  t')  <  t;(t), 

where  v  is  the  valuation  function  defined  as  above.  Thus  in  this  sense  each  remainder  is 
small. 

In  the  boolean  case,  Reference  30,  there  exists  an  n  such  that 

a„ia(tO"  =0 

so  that  the  expression  for  a  becomes 

a  =  Tk  0 

One  useful  application  of  this  result  for  the  boolean  case  is  in  data  compression.  By  encoding 
the  Fj’s  in  run  length  code,  the  image  can  be  represented  by  fewer  bits  of  data,  and  recon¬ 
structed  exactly  once  t  is  known. 

9.  IMAGE  ALGEBRA  AND  MATHEMATICAL  MORPHOLOGY 

This  section  discusses  the  connection  between  the  two  structures  of  the  Image  Algebra 
under  ®  and  Q  and  mathematical  morphology,  in  particular  how  these  structures  general¬ 
ize  morphological  concepts.  Mathematical  morphology,  first  pioneered  by  G.  Matheron, 
Reference  31,  and  J.  Serra,  Reference  32,  depends  on  two  fundamental  set  theoretic  operar 
tions,  namely  the  Minkowski  addition  and  subtraction  of  sets  in  Euclidean  space,  Refer¬ 
ence  33.  For  any  set  ACR“  and  BCR",  Minkowski  addition  and  subtraction  are  defined  as 

A@B  =  {a-l-b :  a€  A,  b€B}  and  ABB  =  (A'eB')'  , 

respectively,  where  B'  =  { -b  :  b  €  B  }  and  A'  denotes  the  complement  of  A  in  R".  The  two 
basic  morphological  operations  of  erosion  and  dilation  of  A  by  B  are  then  defined  as  A©B' 
and  A@B'  respectively.  To  avoid  anomalies  without  practical  interest,  it  is  usually  assumed 
that  B  contains  the  origin.  In  applications  A  represents  a  set  of  pixels  and  B  is  referred  to  as 
a  structuring  element. 
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A  more  general  image  transform,  called  the  Hit  or  Miss  transform,  can  be  defined  in 
terms  of  the  Minkowski  subtraction.  This  transform  is  often  viewed  as  the  universal  mor¬ 
phological  transform  upon  which  the  theory  of  mathematical  morphology  is  based.  Its 
definition  is  as  follows.  Suppose  B  =  (D,E)  is  a  pair  of  structuring  elements,  then  the  Hit  or 
Miss  transform  of  the  set  A  is  given  by  the  expression 

A©B  =  {a:D^CA,  E,,CA'} 

where  denotes  the  translate  D  4-  a  of  D  by  a,  a  £  A,  and  similarly  for  E^.  For  practical 
applications,  it  is  assumed  that  D  H  E  =  <2li.  It  turns  out  that  if  E  =  0,  then  A©B  =  A©D'. 
Thus,  the  Hit  or  Miss  transform  includes  erosion  and,  hence,  dilation  as  a  special  case. 

Of  course,  in  actual  image  processing  we  are  not  dealing  with  arbitrary  sets  in  Euclidean 
space  but  bounded  arrays  of  fixed  size.  Thus,  in  the  above  definitions  we  cannot  simply  sub¬ 
stitute  an  image  for  the  symbol  A,  as  the  expression  A'  would  become  meaningless  to  the 
computer.  What  is  usually  assumed  is  that  A  is  a  set  of  points  within  a  given  image  such  as 
the  set  of  all  black  pixels  in  a  Boolean  image.  Another  observation  is  that  the  Minkowski 
operations  are  not  applied  in  the  same  manner  to  gray  valued  images  as  they  are  to  Boolean 
images.  In  fact,  the  operations  of  erosion  and  dilation  for  gray  valued  images  have  to  be 
redefined  through  the  cumbersome  notion  of  the  umbra  of  an  image,  Reference  32.  In  con¬ 
trast,  the  simple  max  and  min  Image  Algebra  convolutions  defined  by  Q  and  Q  can  be 
used  to  express  the  morphological  operations  of  dilation  and  erosion,  respectively,  for  both 
Boolean  and  gray  valued  images.  In  particular,  if  B  denotes  the  structuring  element,  we  may 
define  a  template  t  =  (T,  <)  corresponding  to  B  by  setting  7(y)  =  B'y,  and  let  ty(x) 
correspond  to  the  values  assigned  to  By'  at  location  x  with  ty(x)  =  0  if  x  T(y).  Then, 
a  □  t  is  equivalent  to  a  dilation  of  A  by  B,  and 
a  □  — t  is  equivalent  to  an  erosion  of  A  by  B,  where  -t  =  (T,-t). 

A  little  care  needs  to  be  taken  in  the  interpretation  of  the  above  two  statements.  For 
instance,  in  the  Image  Algebra  expression  aQ  t,  a  denotes  the  actual  input  image.  Thus,  for 
Boolean  image  processing  (ty(x)  =  0  V  x),  the  output  image  a  Q  t  will  have  the  black  com¬ 
ponents  of  a  (a(x)  =  1)  dilated.  In  contrast,  in  the  morphological  transform  A@B'  A  does 
not  denote  the  input  image  but  only  the  set  of  black  pixels  of  the  input  image,  while  A@B' 
denotes  the  set  of  black  pixels  of  the  output  image.  In  addition,  structuring  elements 
correspond  only  to  invariant  templates  from  X  to  X  with  the  property  that  the  target  pixel 
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is  an  element  of  the  source  configuration  T(y),  a  very  limited  class  of  templates. 

If  t  has  only  zero  weights,  then  t  =  -t  and  a  Q  — t  =  a  Q  t.  Hence  the  eroded  image 
a  Q  — t  is  the  image  obtained  from  a  by  replacing  each  pixel  a(x)  by  the  minimum  value  in 
the  source  configuration  (window)  71[y).  That  this  type  of  min  and  max  operations  on 
Boolean  images  is  equivalent  to  erosions  and  dilations  follows  from  P.  Miller’s  Boolean  alge¬ 
bra  approach  to  morphology,  Reference  30.  Figure  8  represents  the  dilation  aQ  t  of  the 
input  image  a  shown  in  Figure  8(a),  while  Figure  8(c)  represents  the  erosion  of  the  im^e 
a  0 1,  namely,  (a  Q  t)  n  t  which  is  also  called  the  closing  of  a  by  t.  (  Serra,  Reference  32  ). 
The  template  t  used  in  this  example  is  defined  by 


t(x)  = 


The  operational  equivalence  between  0  ,  0  and  dilations  and  erosions  for  gray  valued 

images  stems  from  the  fact  that  when  actually  implementing  dilations  and  erosions  for  image 

processing  purposes,  the  two  equations 

c(x)  =  V(a(x-y)  -I-  b(-y)] 
y 

c(x)=A(a(x-y)-  b(-y)l 
y 

are  used,  where  a(x)  denotes  the  input  pixels  and  b(y)  the  values  of  the  structuring  element 
B  (Serra,  pg.  441,  Reference  32  ). 

In  the  Boolean  case  it  is  a  simple  exercise  to  show  that  the  Hit  or  Miss  transform  can  be 
expressed  as 

A©B  =(A©D')n(A'©E'). 

Defining  t  =  ( T,  t)  by  7(y)  =  D^',  ty(x)  =  0  and  s  =  ( 5.  s)  by  5(y)  =  E^'  and  Sy(x)  =  0.  we 
obtain  the  equivalent  Image  Algebra  expression 


However,  there  is  an  even  simpler  Image  Algebra  formulation  of  the  Hit  or  Miss 
transform  which  does  not  employ  the  notions  of  minimum  or  erosion.  Suppose  T(y)  = 

{xi,  ....  xj  and  5(y)  =  {x^+i,  .  .  .  ,x„}.  Then 

k 

Xm(a@r),  where  m  = 

is  equivalent  to  computing  the  Hit  or  Miss  transform,  where  r  =  (.R,  r)  is  defined  by  J2(y)  = 
T(y)  U  5(y)  and  ry(xi)  =  2^"^ 

Since  Image  Algebra  can  express  the  two  basic  morphological  operations  of  erosion  and 
dilation,  it  is  clear  that  the  subalgebra  A  =  (R^,  Rx-»x>  +-  *>  V,  Q  )  of  the  full  Image  Alge¬ 
bra  includes  mathematical  morphology  as  a  special  case.  However,  even  the  subalgebra  A  is 
a  much  larger  algebra  than  that  provided  by  the  algebra  defined  by  mathematical  morphol¬ 
ogy.  Templates  are  more  general  objects  than  structuring  elements.  They  can  vary  in  size, 
shape  and  weights  from  point  to  point  and  they  can  map  between  different  structures  and 
dimensions  if  we  replace  R  by  F  and  Rx-,x  ^Y-»x-  Thus,  an  expression  of  form  aQ  t 
can  express  a  far  more  complex  process  than  a  simple  dilation.  For  example,  if  a  denotes  the 
input  image  shown  in  Figure  36  and  t  =(T,  t)  is  defined  by  7(i,j)  =  {(x,y)}  with 


1 


lo 


2  2 

if  0.9  <  ~  +  <1.1 

p-  q2 

or  0.9  <  -  -2^  <  1.1 

d^  q" 


otherwise 


where  q=15,  p^  =  q"  +  c^,  d^  =  c^  —  q‘  and  c  =  30,  then  a  0 1  is  obviously  not  a  dilation  nor 
is  it  expressible  in  terms  of  dilations  and  erosions  when  starting  with  the  input  image  a.  The 
input  and  output  images  are  shown  in  the  next  two  figures. 
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Figure  36.  Input  Image  a 


Figure  37.  Image  a  Q  t 


In  the  Boolean  case,  morphological  operators  can  be  accomplished  with  either  ®  or  Q 
by  observing  that  for  a  Boolean  image  a  and  a  template  t  with  ty(x)  =  0  V  x  €  T(y), 
a  El  t  =  a®  t,  where  t  —  (T,  f)  and  i^(x)  =  1  if  x  €  T{y). 
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SECTION  V 

IMAGE  ALGEBRA  SOFTWARE  DEVELOPMENT 


The  goal  of  the  Image  Processing  Language  Project  is  the  development  of  an  algebra 
suitable  for  describing  image  processing  transformations.  A  successful  Image  Algebra  can 
serve  as  a  direct  representation  of  image  processing  algorithms  if  embedded  into  a  program¬ 
ming  language.  Since  the  early  part  of  Phase  I  of  the  contract,  we  have  been  conducting 
work  in  the  implementation  of  the  Image  Algebra  in  FORTRAN.  In  this  section,  we  give  a 
brief  chronology  of  that  work,  as  well  as  other  work  in  software  development  that  we  have 
carried  out  during  the  duration  of  this  contract. 

A  preprocessor  translates  a  source  program  in  an  extended  version  of  some  language  L 
into  a  program  in  L.  Since  many  image  processing  application  programs  are  written  in 
FORTRAN,  an  Image  Algebra  FORTRAN  (lAF)  preprocessor  was  undertaken  in  the  Spring 
of  1984  as  a  senior  project  by  W.K.  Perry.  The  LAF  preprocessor  translates  FORTRAN 
programs  with  image  algebra  extensions  into  standard  FORTRAN  programs.  This  gave 
users  the  ability  to  define  translation  invariant  templates  and  to  execute  binary  image-image 
and  image-template  operations  using  a  notation  similar  to  the  Image  Algebra.  Various 
ASCII  characters  (such  as  I ,  # ,  .  etc.)  were  used  in  this  original  version  of  Image  Alge¬ 

bra  FORTRAN  to  denote  the  Image  Algebra  symbols  like  ® ,  and  |2  ,  that  do  not  appear 
on  a  standard  keyboard. 

During  the  Fall  of  1985,  work  on  the  preprocessor  led  to  the  development  of  a  notation 
for  defining  translation  variant  templates.  This  required  a  substantial  change  to  the  design 
of  the  preprocessor.  A  new  program  unit,  the  template  definition,  was  added  to  the  Image 
Algebra  Fortran  language. 

Until  the  Fall  of  1985,  we  had  been  using  a  Printronix  line  printer  as  our  primary  image 
output  device.  Our  gray-toning  software  permitted  us  to  print  images  on  the  dot  matrix 
printer  using  six  by  six  squares  of  dots  to  represent  a  single  pixel  in  an  image.  At  this  time, 
our  department  acquired  a  300  dot  per  inch  resolution  laser  printer  and  we  modified  our 
gray-toning  software  to  use  the  laser  printer  as  an  output  device.  This  improved  the  resolu¬ 
tion  of  our  image  output  by  a  factor  of  about  16. 

During  this  period,  the  LAF  preprocessor  came  into  heavy  use  since  the  preprocessor  had 
been  distributed  to  a  number  of  companies  who  had  expressed  an  interest  in  having  this 
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Image  Algebra  software  tool.  It  became  clear  that  although  the  preprocessor  was  adequate 
to  the  expression  of  both  Image  Algebra  expressions  and  much  of  standard  FORTRAN,  it 
was  not  ideal  at  combining  both  these  tasks  within  a  single  program.  A  number  of  shortcom¬ 
ings  were  identified: 

(1)  Image  Algebra  expressions  could  only  appear  within  the  context  of  a  FORTRAN 
assignment  statement, 

(2)  although  arbitrarily  many  image-image  operations  using  FORTRAN  arithmetic 
operations  could  appear  in  a  single  assignment  statement,  if  an  image-template  operation  or 
non-FORTRAN  image-image  operation  appeared  in  an  assignment  statement,  it  must  be  the 
only  operation  in  that  assignment  statement, 

(3)  the  set  of  variable  names  used  by  the  preprocessor  could  clash  with  many  user 
defined  names, 

(4)  although  the  image  expression  of  the  form  A  +  1  would  be  acceptable,  1  +  A 
would  not  be  recognized  as  an  image  expression, 

(5)  an  image  could  appear  both  on  the  left  hand  side  of  an  assignment  statement  and  as 
an  operand  of  an  image-template  operation  on  the  right  hand  side, 

(6)  no  statement  involving  declaration  of  arrays  or  assignments  to  images  could  be  con¬ 
tinued  on  multiple  lines,  and 

(7)  almost  no  error  checking  was  performed  by  the  preprocessor. 

These  problems  stemmed  from  the  fact  that  the  preprocessor  performed  a  simple  kind  of 
pattern  matching  to  drive  its  translation.  This  simple  pattern  matching  method  is  much  too 
weak  to  handle  the  subtleties  of  the  FORTRAN  syntax.  Although  various  improvements 
were  made  to  the  preprocessor,  it  retained  this  simple  basis,  making  the  improvement  of 
many  of  its  deficiencies  impossible  without  a  complete  redesign  effort. 

By  the  Spring  of  1986,  it  had  become  clear  to  us  that  the  Image  Algebra  is  readily 
implementable  in  a  programming  language  and  can  yield  fantastic  benefits  in  the  speeding  up 
of  the  program  development  process.  We  had  demonstrated  the  lAF  preprocessor  several 
times  at  program  reviews  and  met  with  enthusiasm  and  great  interest  on  the  part  of  practi¬ 
cally  all  of  those  who  had  seen  its  use  in  the  development  of  code.  In  preparation  for  poten¬ 
tial  military  funded  development  of  the  Image  Algebra  in  a  production  quality  setting,  we 
contacted  the  several  Ada  compiler  producers.  The  companies  contacted  produced  Ada 
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compilers  for  a  wide  variety  of  machines  at  that  time,  their  products  had  good  reputations. 
We  felt  that  it  was  important  to  contact  Ada  compiler  producers  to  determine  the  feasibility 
of  cooperative  development  of  an  Image  Algebra  Ada  compiler,  that  is,  a  compiler  for  Ada 
that  provided  the  user  with  Image  Algebra  operations  and  operands.  Such  a  cooperative 
effort  with  a  company  having  existing  Ada  expertise  could  yield  a  workable  Ada  product 
with  less  lead  time  than  development  from  the  ground  up.  While  all  the  companies  con¬ 
tacted  indicated  that  cooperative  work  could  be  carried  out,  their  estimates  of  price  and  level 
of  interaction  differed.  Estimates  of  cost  of  release  of  source  code  varied  from  $150,000  to 
$1,000,000. 

During  the  summer  of  1986,  we  were  able  to  arrange  a  loan  of  a  Sun  3/160C  worksta¬ 
tion.  At  that  time,  we  were  able  to  develop  an  image  display  package  callable  by  lAF  pro¬ 
grams.  This  gave  us  the  capability  of  displaying  pseudo-color  representations  of  8-bit 
imagery  during  processing.  We  obtained  our  own  Sun  Workstations  via  a  DoD  University 
Research  Initiative  Program  grant  in  December  of  1986  and  were  then  able  to  use  the 
preprocessor  to  develop  nontrivial  example  programs  since  our  throughput  time  was  greatly 
enhanced  by  the  workstation’s  real-time  image  display  capability. 

While  the  Image  Algebra  FORTRAN  preprocessor  provides  us  with  the  ability  to  test 
concepts  of  the  Image  Algebra,  it  does  not  provide  a  rich  environment  for  the  testing  of 
Image  Algebra  implementation  concepts.  A  preprocessor  has  inherent  limitations  that  affect 
its  usability. 

•  Runtime  efficiency  is  poor  since  the  language  extensions  provided  by  a  preprocessor 
must  be  implemented  in  terms  of  a  target  high-level  language  that  was  deficient  in  pro¬ 
viding  the  desired  extensions  in  the  first  place.  Hence,  the  extensions  may  not  be  imple¬ 
mented  in  the  most  appropriate  way  on  any  given  machine. 

•  Translation  efficiency  in  a  preprocessor  is  poor  since  an  extra  translation  step — 
extended  host  language  to  host  language — must  be  carried  out  to  produce  a  running 
program. 

•  User  comprehension  is  degraded  by  having  to  understand  one  extra  interface  level. 
Runtime  error  messages  are  often  reported  with  respect  to  the  translated  host  language 
version  of  the  program  rather  than  the  source  language. 

Bearing  in  mind  these  limitations  of  the  preprocessor  approach  to  language  implementar 
tion,  in  the  fall  of  1986,  we  decided  to  study  the  issues  involved  in  implementing  the  Image 
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Algebra  in  a  compiler.  A  compiler  is  a  program  that  translates  a  high  level  language  source 
program  into  an  executable  machine  language  object  program.  Such  a  language  implementa¬ 
tion  does  not  suffer  from  the  problems  cited  above  for  preprocessors.  A  typical  compiler  con¬ 
sists  of  two  phases:  the  front  end,  which  is  a  machine  independent  compilation  phase 
translating  the  source  code  of  the  compiled  language  into  an  intermediate  code;  and  the  back 
end,  which  is  a  machine  dependent  phase  translating  the  intermediate  code  into  the  machine 
language  of  the  target  machine. 

One  question  we  addressed  at  this  time  was  what  kind  of  intermediate  code  might  be 
appropriate  for  use  in  such  a  compiler.  After  a  review  of  the  relevant  compiler  literature  we 
decided  that  since  Stanford  U-Code  haul  been  used  in  the  development  of  a  machine  indepen¬ 
dent  global  data  flow  optimizer  34  it  might  be  a  suitable  intermediate  code  for  the  Image 
Algebra.  Global  data  flow  analysis  is  the  study  of  the  flow  of  data  values  through  programs. 
Many  sorts  of  code  improvement  techniques  can  be  performed  using  information  derived 
from  the  analysis  of  the  dataflow  of  a  program.  Since  execution  of  a  single  Image  Algebra 
operation  typically  involves  many  thousands  or  even  millions  of  operations  on  data  in  a 
highly  structured  way,  there  are  many  opportunities  to  exploit  global  dataflow  properties 
during  the  translation  of  source  code  to  machine  executable  object  code. 

During  this  period  W.K.  Perry  undertook  initial  steps  required  in  the  development  of  a 
compiler  for  the  C  language- with  Image  Algebra  extensions  35.  The  designed  compiler  con¬ 
sists  of  three  major  pieces  of  code:  a  front  end,  translating  Image  Algebra  C  into  assembly 
code;  an  assembler,  translating  the  assembly  code  into  object  code;  and  a  machine  simulator 
interpreting  the  produced  object  code.  We  avoided  compiling  directly  to  the  machine 
language  of  any  particular  machine  so  that  the  implementation  could  be  ported  to  machines 
other  than  the  one  on  which  it  was  originally  implemented.  This  work  demonstrated  that 
U-Code,  designed  for  a  uniprocessor  model  of  computation,  is  in  fact  a  weak  language  for 
implementing  Image  Algebra  operations.  A  vector  variant  of  U-Code,  V-Code,  was 
developed  to  permit  more  efficient  execution  of  operations  in  the  compiler-assembler- 
simulator  environment. 

Another  student  began  a  study  of  global  dataflow  based  code  improvements  suggested 
by  the  kind  of  code  produced  by  straightforward  implementation  of  the  Image  Algebra  in  a 
uniprocessor  setting  36.  Her  initial  task  was  to  study  existing  Image  Algebra  algorithms; 
determine  what  sorts  of  data  interdependencies  appeared  in  the  generated  code;  and 
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determine  techniques  of  code  modification  that  could  be  mechanically  applied  to  yield  code 
with  identical  results  but  faster  execution  speed.  One  of  the  initial  algorithms  she  chose  to 
look  at  was  histogram  computation.  Suppose  a  is  an  image.  One  algorithm  for  computing 
the  histogram  h  of  a  using  Image  Algebra  operations  is  the  following: 

for  i  in  Aa  to  Va  do 

h(.)  -  Eu(»)) 


Expressed  in  Image  Algebra  FORTRAN,  we  have  the  following  algorithm: 

do  10  Is  !mln  a,  'max  a 

h(l)  =  !sum  (  a  =  1  ) 

10  continue 

This  leads  to  generation  of  the  following  FORTRAN  code  in  the  current  preprocessor: 
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*  do  10  1=  Imln  a,  Imax  a 
XXSI01=A(1,1) 

DO  99000  XXLC:V2=1.NCQLS 
DO  99001  XXLCV1=1,NR0WS 
IF  (A(XXLCV1.XXLCV2)  LT.XXSIOl)  THEN 
XXSI01=A(XXLCV1 , XXLCV2) 

ENDIF 

99001  CONTINUE 
99000  CONTINUE 

XXSI02=A(1.1) 

DO  99002  XXLCV2=l,NCaLS 
DO  99003  XXLCV1=1,NR0WS 
IF  (A (XXLCVl , XXLCV2) . GT . XXSI02)  THEN 
XXSI02=A (XXLCVl , XXLCV2) 

ENDIF 

99003  CONTINUE 

99002  CONTINUE 

DO  10  I=XXSI01,XXSI02 

*  h(l)  =  !sum  (  a  =  i  ) 
XXIND3=1 

DO  99004  XXLCV2=-1,NCQLS 

DO  99005  XXLCV1=.' .NRQWS 

IF  (A (XXLCVl .  yXLCJ2) . EQ . I)  THEN 

XXI001(XXIND3)=1 

ELSE 

XXIOOl  (XXIND3)=0 
ENDIF 

XXIND3=XXIND3+1 

99005  CONTINUE 

99004  CONTINUE 
XXSI03=0 
XXIND2=1 

DO  99006  XXLCV2=1.NCQLS 
DO  99007  XXLCV1=1, NRQWS 
XXSI03=XXSI03+XXI001 (XXIND2) 
XXIND2=XXIND2+1 
99007  CONTINUE 

99006  CONTINUE 
H(I)=XXSI03 

10  continue 
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One  of  the  optimizations  being  investigated  is  that  of  loop  joining.  Loop  joining  consists 
roughly  of  concatenating  the  bodies  of  loops  together  if  they  have  the  same  loop  bounds. 

This  cam  only  be  done  if  certain  properties  of  the  loops  hold  that  guarantee  that  the  joining 
will  not  change  the  behavior  of  each  of  the  loops.  The  code  as  modified  by  loop  joining  is  the 
following; 

*  do  10  1=  Imln  a.  Imaix  a 
XXSI01=A(1.1) 

C  Assignment  to  XXSI02  hoisted  up  from  second  loop 
XXSI02=A(1 . 1) 

DO  99000  XXLCV2=l.NCaLS 
DO  99001  XXLCV1=1,NR0WS 
IF  (A(XXLCV1,XXLCV2) .LT.XXSIOl)  THEN 
XXSI01=A(XXLCV1 ,XXLCV2) 

END  IF 

C  Code  that  follows  was  body  of  second  loop 

IF  (ACXXLCV1,XXLCV2) .GT.XXSI02)  THEN 
XXSI02=A (XXLCVl , XXLCV2) 

END  IF 

99001  CONTINUE 
99000  CONTINUE 

DO  10  I=XXSI01,XXSI02 

♦  h(l)  =  !s\im  (  a  =  1  ) 

XXIND3=1 

C  These  two  assignments  hoisted  up  from  second  loop 
XXSI03=0 
XXIND2=1 

DO  99004  XXLCV2=1,NC0LS 

DO  9900S  XXLCV1=1.NR0WS 

IF  (A(XXLCV1,XXLCV2) .EQ.I)  THEN 

XXI001CXXIND3)=1 

ELSE 

XXI001(XXIND3)=0 
END  IF 

XXIND3=XXIND3+1 

C  code  that  follows  was  body  of  the  second  loop 
XXSI03=XXSI03+XXI001 <XXIND2) 

XXIND2=XXIND2+1 
99007  CONTINUE 
99006  CONTINUE 

HCly  =XXSI03 
10  continue 
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Noting  that  the  two  induction  variables  XXIND3  and  XXIND2  have  the  same  value,  we 
can  perform  what  we  call  backward  code  motion  to  yield  the  following  program: 

*  do  10  1=  !mln  a,  !max  a 
XXSI01=A(1.1) 

XXSI02=AC1.1) 

DO  99000  XXLCV2=l.NCaLS 
DO  99001  XXLCV1=1,NR0WS 
IF  (A(XXLCV1.XXLCV2) .LT.XXSIOl)  THEN 
XXSI01=A (XXLCVl . XXLCV2) 

END  IF 

IF  (A (XXLCVl , XXLCV2) . GT . XXSI02)  THEN 
XXSI02=A (XXLCVl .XXLCV2) 

END  IF 

99001  CONTINUE 
99000  CONTINUE 

DO  10  I=XXSI01.XXSI02 

*  h(l)  =  !s\im  (  a  =  1  ) 

XXIND3=1 

XXSI03=0 

DO  99004  XXLCV2=1,NC0LS 
DO  99005  XXLCV1=1.NR0WS 
IF  (A  (XXLCVl .  XXLCV2)  .  EQ .  I)  THEN 
XXI001(XXIND3)=1 

C  assignment  to  XXSI03  was  moved  back  into  this  if  branch 
XXSI03=XXSI03+XXI001 (XXIND3) 

ELSE 

XXI001(XXIND3)=0 

C  Assignment  to  XXSI03  was  moved  back  into  this  if  branch 
XXSI03=XXSI03+XXI001 (XXIND3) 

ENDIF 

XXIND3=XXIND3+1 
99007  CONTINUE 
99006  CONTINUE 

H(I)=XXSI03 
10  continue 
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Copy  propagation  is  the  removing  of  assignments  to  variables  whose  values  are  copied 
via  assignment  to  other  variables  and  never  used  again.  Performing  copy  propagation,  we 
are  left  with  this  program: 

*  do  10  1=  Imln  a.  !max  a 
XXSI01=A(1.1) 

XXSI02=A(1.1) 

DO  99000  XXLC:V2=1.NC0LS 
DO  99001  XXLCV1=1.NR0WS 
IF  (A00CLCV1.XXLCV2)  .LT.XXSIOl)  THEN 
XXSI01=A(XXLCV1 . XXLCV2) 

ENDIF 

IF  (A(XXLC:V1.XXLC:V2)  .GT.XXSI02)  THEN 
XXSI02=A (XXLCVl , XXLCV2) 

ENDIF 

99001  CONTINUE 
99000  CONTINUE 

DO  10  I=XXSI01,XXSI02 

*  hCl)  =  !s\im  (  a  =  1  ) 

C  XXIND3  has  been  removed  because  it  now  has  no  effect 
C  on  the  behavior  of  the  program 

C  XXSI03  has  been  removed  since  it  is  used  only  to  directly  assign 
C  it’s  value  to  H(I) 

H(I)  =  0 

DO  99004  XXLCV2=1.NC0LS 
DO  9900S  XXLCV1=1.NR0WS 
IF  (A (XXLCVl . XXLCV2) . EQ . I)  THEN 

C  XXIOOl  has  been  removed  since  it  is  used  only  to  directly  assign 
C  it’s  value  to  H<I)  . 

H(I)=HCI)+1 

ELSE 

H(I)  =  H(I)+0 
ENDIF 

99007  CONTINUE 
99006  CONTINUE 
10  continue 
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We  can  now  modify  by  simplification  of  algebraic  expressions.  This  makes  the  assign¬ 
ment  of  H(I)  +  0  to  H(I)  a  needless  operation,  thus  we  can  remove  the  else  branch  con¬ 
taining  this  statement,  yielding  the  following  code: 

*  do  10  1=  !mln  a,  'max  a 
XXSI01=A(1,1) 

XXSI02=A(1,1) 

DO  99000  XXLCV2=1,NC0LS 
DO  99001  XXLCV1=1.NR0WS 
IF  CA(XXLCV1.XXLCV2) .LT.XXSIOl)  THEN 
XXSI01=A(XXLCV1 , XXLCV2) 

END  IF 

IF  (A(XXLCV1,XXLCV2) .GT.XXSI02)  THEN 
XXSI02=A (XXLCVl . XXLCV2) 

ENDIF 

99001  CONTINUE 
99000  CONTINUE 

DO  10  I=XXSI01.XXSI02 

*  h(i)  =  Isiim  C  a  =  1  ) 

HCI)  =  0 

DO  99004  XXLCV2=l,NCaLS 
DO  99005  XXLCV1=1,NRQWS 
IF  (A (XXLCVl . XXLCV2) . EQ . I)  THEN 
H(I)=HCI)+1 

C  Else  branch  with  no  effect  removed 
ENDIF 

99007  CONTINUE 
99006  CONTINUE 
10  continue 
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Loop  Splitting  consists  of  breaking  the  statements  of  one  loop  apart  and  placing  them 
into  two  loops.  This  can  be  useful  if  it  is  possible  to  further  optimize  one  of  the  loops  in  the 
absence  of  the  statements  in  the  other.  Following  this  idea,  we  split  the  loop  over  the  gray 
value  range  of  the  image  a,  yielding  this  version  of  the  program: 

*  do  10  1=  Imln  a,  Imax  a 
XXSI01=A(1.1) 

XXSI02=ACl . 1) 

DO  99000  XXLC:V2=1.NC0LS 
DO  99001  XXLCV1=1,NRQWS 
IF  (A(XXLCV1.XXLCV2) .LT.XXSIOl)  THEN 
XXSI01=A (XXLCVl , XXLCV2) 

ENDIF 

IF  (A (XXLCVl . XXLCV2) . GT . XXSI02)  THEN 
XXSI02=A (XXLCVl , XXLCV2) 

ENDIF 

99001  CONTINUE 
99000  CONTINUE 

C  This  loop  has  been  introduced 
DO  5  I=XXSI01,XXSI02 
H(I)  =  0 
5  CONTINUE 

DO  10  I=XXSI01.XXSI02 

*  h(l)  =  !sum  (  a  ==  1  ) 

DO  99004  XXLCV2=1.NC0LS 

DO  99005  XXLCV1=1.NR0WS 

IF  (A (XXLCVl , XXLCV2) . EQ . I)  THEN 

H(I)=H(I)+1 

ENDIF 

99007  CONTINUE 
99006  CONTINUE 
10  continue 
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Given  certain  conditions  on  the  body  of  a  nested  loop,  it  is  possible  to  interchange  the 
loop  iterations.  In  this  version  of  the  program,  we  interchange  the  loop  over  I  with  the  two 
loops  nested  within  it; 

*  do  10  1=  !mln  a.  !max  a 
XXSI01=A(1.1) 

XXSI02=A(1 . 1) 

DO  99000  XXLCV2=1.NCQLS 
DO  99001  XXLCV1=1.NR0WS 
IF  CA(XXLCV1.XXLCV2) .LT.XXSIOl)  THEN 
XXSI01=A (XXLCVl . XXLCV2) 

END  IF 

IF  (A  OOCLCVl .  XXLCV2)  .  GT .  XXSI02)  THEN 
XXSI02=A  (XXLCVl .  XXLC::V2) 

ENDIF 

99001  CONTINUE 
99000  CONTINUE 

DO  S  I=XXSI01.XXSI02 
H(I)  =  0 
S  CONTINUE 

*  h(l)  =  Isxim  (  a  =  1  ) 

DO  99004  XXLCV2=1 . NCOLS 

DO  9900S  XXLCV1=1.NR0WS 

C  loop  on  I  has  been  shifted  to  be  innermost 
DO  10  I=XXSI01.XXSI02 
IF  (A(XXLCV1.XXLCV2) .EQ.I)  THEN 
H(I)=H(I)+1 
ENDIF 

10  continue 
99007  CONTINUE 
99006  CONTINUE 
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If  two  consecutive  if  statements  without  else  branches  are  mutually  exclusive,  that  is,  if 
given  that  one  of  the  if  tests  is  satisfied  the  other  is  not  satisfied  and  vice  versa,  then  we  can 
join  these  statements  into  an  if-then-else  statement.  Likewise,  if  the  body  of  a  do-loop  con¬ 
sists  of  only  a  single  if-statement  with  no  else  branch,  and  the  if  condition  is  satisfied  only 
when  the  loop  index  takes  on  a  particular  value  v,  then  the  loop  can  be  removed  and  all 
references  to  the  do-loop  index  within  the  if  can  be  replaced  by  v.  Performing  these 
modifications  results  in  the  following  program; 

*  do  10  1=  !mln  a,  !max  a 
XXSI01=A(1.1) 

XXSI02=A(1.1) 

DO  99000  XXLCV2=1.NC0LS 
DO  99001  XXLCV1=1.NR0WS 
IF  <A(XXLC:V1,XXLCV2)  .LT.XXSIOl)  THEN 
XXSI01=A(XXLCV1 ,XXLCV2) 

C  We  join  this  IF  with  the  previous  one  by  else  introduction 
ELSE  IF  (A(XXLCV1,XXLCV2) .GT.XXSI02)  THEN 
XXSI02=A (XXLCVl . XXLCV2) 

ENDIF 

99001  CONTINUE 

99000  CONTINUE 

DO  5  I=XXSI01,XXSI02 
HCI)  =  0 
5  CONTINUE 

*  hCl)  !sum  C  a  =  1  ) 

DO  99004  XXLCV2=1,NC0LS 

DO  99005  XXLCV1=1.NR0WS 

C  This  statement  was  guarded  by  a  loop  joinable  IF 

H (A (XXLCVl , XXLCV2) ) =H (A (XXLCVl . XXLCV2) ) + 1 

99007  CONTINUE 

99006  CONTINUE 

At  this  point,  the  example  has  been  modified  to  the  point  where  it  is  effectively  no 
different  from  the  sort  of  algorithm  that  a  programmer  might  develop  by  hand  for  histogram 
computation.  The  speedup  achieved  by  this  modification  is  about  a  factor  of  300  on  our  Sun 
Microsystems  3/280  uniprocessor.  We  are  continuing  our  study  of  code  improvement  tech¬ 
niques,  trying  to  identify  necessary  and  sufficient  conditions  to  insure  that  the  modified  code 
will  run  faster  than  the  original  code.  We  are  also  searching  for  ways  to  decrease  the 
amount  of  time  needed  to  determine  the  best  restructuring  of  the  original  code. 
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It  is  our  hope  that  this  seminal  work  on  implementation  of  the  Image  Algebra  can  serve 
as  a  starting  point  for  the  development  of  production  quality  software  products  providing 
computer  vision  program  developers  with  the  capability  to  solve  image  processing  tasks  with 
directly  implemented  Image  Algebra  primitives.  Our  experience  has  shown  that  not  only  is 
it  possible  to  enhance  the  productivity  of  programmers  using  the  Image  Algebra,  but  it  is 
also  possible  to  efficiently  implement  Image  Algebra  primitives  in  a  programming  language 
system. 
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